

OBSERVATIONAL AND TH EOR ETICAL INVESTIGATIONS 

JN SOLAR SEISMOLOGY „ „ ^ 

- m- 


NASA Grant NAG W- 1677 


Semiannual Reports No. 1, 2, 3, and Final Report 

For the peri od 1 March 1 giQhro'ugh 3T August 1 990 


0 ^/ 7 ^ 

p- v* 


For the period 1 March 1989 thro ugh 28 February 1991 

Principal Investigator 


Robert W. Noyes 



March 1 992 


Prepare d for 



National Aeronaut i c s and Space Administration 


« - Washington, D.C. 20546 


Smithsoni an Institution 
Astrophysical Observatory 
Cambridge, M assachu setts 02138 


The Smithsonian Astro physical Observatory 
is a member of the 

Harvard-Smithsonian Center for Astrophysics 



The NASA Technical Officer for this Gra nt is Dr. William J. Wagner, Solar Physics 
Branch, Space Physics Division, Code SS, NASA Headquarters, Washington, D.C. 20546 


(NASA-CR -190245) OBSERVATIONAL AND 
THEORETICAL INVESTIGATIONS IN SOLAR 
SEISMOLOGY Final Semiannual Reports No. 1-3, 
1 Mar. 1989 - 31 Aug. 1990 and 1 Mar. 1989 - 
23 Feb. 1991 (Smithsonian As trophys i cal 


N 9 2-26247 


Unc 1 as 

G 3/92 0085174 


J 




Observational and Theoretical Investigations 
in Solar Seismology 

Grant Number NAGW-1677 


I. INTRODUCTION 

This is the final report on a project to develop a theoretical basis for interpreting solar 
oscillation data in terms of the interior dynamics and structure of the Sun. Work funded 
under this project was carried out principally by Dr. Michael Ritzwoller, Dr. Pawan 
Kumar, and Dr. Sylvain Korzennik, under the overall guidance of the Principal Investiga- 
tor, from March 1989 through August 1991. In sections D, in, and IV we discus: heir 
three related areas of work, and in section V give references to four papers written with 
support from this grant. An Appendix gives complete copies of two of these papers and 
abstracts and summaries of the other two. 

n. STUDIES OF THE HELIOSEISMIC SIGNATURES OF DIFFERENTIAL ROTA- 
TION AND CONVECTION IN THE SOLAR INTERIOR. 

The now-traditional method to infer differential rotation in the solar interior is to expand 
the rotation rate Q(r,0) as a power series in cos k 0 or in P k (cos0), and relate the expansion 
coefficients to the observed helioseismic mode splittings. These splittings may themselves 
may be expanded in Legendre polynomials of m//, where m and / are the azimuthal order 
and degree, respectively, of the oscillation modes. This expansion has the appeal that 
because different terms of the expansion are orthogonal over the sphere, observational 
errors in determining the different terms of the expansion are nearly independent. The 
internal angular velocity, or differential rotation, is then determined by inversion of the 
splitting coefficients. This approach has the difficulties that (1) it is difficult to generalize 
to non-axisymmetric flows, (2) it is computationally cumbersome, and (3) there is cross- 
talk between the various terms in the expansion of Q. 

Ritzwaller and Eugene Lavely, with partial support from this grant, have developed a uni- 
fied approach to the helioseismic forward and inverse problem of differential rotation 
(Ritzwoller and Lavely 1991). In this approach the differential rotation is represented as 
the axisymmetric component of a more general toroidal flow field. A better choice of 
basis functions for differential rotation allows determination of a set of vector spherical 
harmonic expansion coefficients for the rotation that are decoupled so that each degree of 
differential rotation can be estimated independently from all other degrees. 

Lavely and Ritzwoller, also with partial support from this grant, have carried out a funda- 
mental study of the effect of global-scale steady-state convection on helioseismic oscilla- 
tions (Lavely and Ritzwaller 1992). They have derived the basic theory governing the 
influence of convection and associated structural asphericities on oscillation frequencies, 
without the usual assumptions of an axisymmetric model. They represention the eigen- 
functions of a reference spherical model with vector spherical harmonics, and employed 
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quasi-degenrate perturbation theory to derive general matrix elements governing mode 
coupling and splitting caused by convection and structural asphericities. This formalism 
may be applied to models of giant-cell convection, allowing determination of whether such 
flows bias recently estimated differential rotation profiles. Or, it could be used to test 
hypothesized pole-equator differences in temperature near the top of the convection zone 
as a source of observed even-degree splitting coefficients. 

A complete copy of the first of the two above-referenced papers, and the abstract and 
summary of the second, are included in the appendix to this Final Report. 

m. WAVE GENERATION BY TURBULENT CONVECTION. 

Dr. Pawan Kumar, supported in part by this grant, has studied the generation of solar p- 
mode oscillations by turbulent convection in collaboration with P. Goldreich (Goldreich 
and Kumar, 1990). They used a simplified model of an adiabatically stratified upper solar 
convection zone overlaid by a convectively stable isothermal atmosphere, and studied the 
rate at which convective energy is converted into energy of trapped p-modes, f-modes, g- 
modes, and travelling acoustic waves. They found that wave generation is concentrated at 
the top of the convection zone where the turbulent Mach number M t peaks. The effi- 
ciency tj of power input into trapped p and f modes, and into travelling acoustic waves, 
was found to vary as 'n~M t 7 ' 5 , and that into g modes was found to vary as M, For p- 
modes the energy input depends on frequency as co a where a=(2m 2 +7 m-3)/ (m+3) , and m 
is the polytropic index. This agrees with the observed finding of Libbrecht (1988,ESA 
SP-286, p 3) that the power input varies as to 8 , if we set m=4; in fact this value of m is 
close to the polytropic index that fits the density profile in the upper solar convection 
zone. The agreement is strong evidence that wave emission by turbulent convection actu- 
ally is the process by which p-modes are excited in the solar atmosphere. 

Also, Dr. Kumar, with partial support by this grant, collaborated in preparation of a review 
of theories of excitation of oscillation modes in the Sun (Cox et al 1991). This review 
contains a summary of Kumar and Goldreich’ s work on 3-mode coupling, which shows 
that such mode couplings, which had previously been considered a good candidate for lim- 
iting the energies of overstable p-modes, are not strong enough for that purpose. That 
result casts doubt on the idea that the modes are excited by overstability, and therefore 
provides further support for mode excitation by coupling with the motions of convection. 

Copies or summaries of these two papers are also included in the Appendix. 

IV. STUDY OF ANTIPODAL SUNSPOT IMAGING AN ACTIVE REGION TOMOG- 
RAPHY 

During the last stages of this grant, we were joined by Dr. Sylvain Korzennik. Dr Korzen- 
nik has investigated, under partial support of the grant, the signatures of magnetic field 
structures in helioseismic data. Specifically, he carried out a search in Ml Wilson data 
for the helioseismic signature of active regions and sunspots at their antipodal point; 
analogous signals have been seen in geo-seismology, and might be expected in 
helioseismology, because (a) acoustic power is known to be absorbed in sunspots, and (b) 
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helioseismic waves emanating from a point should interfere constructively at its antipodal 
position if the wave lifetime is long enough to establish global coherence and the 
wavelength short enough to prevent blurring of the antipodal "image". Korzennik’s 
analysis of high wavenumber data (2 arcsec per pixel) showed no observable antipodal 
signature from several data sets where moderate-size sunspots were known to exist on the 
invisible hemisphere of the sun. This suggests that there is a lack of global coherence of 
high wavenumber acoustic modes. However, the analysis is still in process. 

The observed acoustic energy deficit in active regions suggests that "active region tomog- 
raphy" will be an important tool for study of the subsurface structure of magnetic active 
regions, particularly with very high-resolution (high wavenumber) data such as may be 
expected from the SOI investigation on SOHO. Dr. Korzennik, with support from this 
grant, began an effort to develop the theoretical and analytical tools to characterize the 
acoustic field fully, with particular attention to the phase relations between velocity and 
intensity perturbations associated with the waves. This work is now being continued, with 
support from the NASA SOI investigation (P. Scherrer, PI), of which Professor Noyes is a 
co-investigator. 

Publication of the results of Dr. Korzennik’s work on antipodal imaging of sunspots and 
active region tomography is currently in preparation. 
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ABSTRACT 

We present a general, degenerate perturbation theoretic treatment of the helioseismic forward and inverse 
problems for solar differential rotation. Our approach differs from previous work in two principal ways. First, 
in the forward problem, we represent differential rotation as the axisymmetric component of a general toroidal 
flow field using vector spherical harmonics. The choice of these basis functions for differential rotation over 
previously chosen ad hoc basis functions (e.g., trigonometric functions or Legendre functions) allows the solu- 
tion to the forward problem to be written in an exceedingly simple form leqs. [32]-[37]). More significantly, 
their use results in inverse problems for the set of radially dependent vector spherical harmonic expansion 
coefficients, which represent rotational velocity, that decouple so that each degree of differential rotation can 
be estimated independently from all other degrees (eqs. [56] & [61]-[63]). Second, for use in the inverse 
problem, we express the splitting caused by differential rotation as an expansion in a set of orthononnal poly- 
nomials that are intimately related to the solution of the forward problem (eqs. [5] and [54]). The orthonor- 
mai polynomials are Clebsch-Gordon coefficients and the estimated expansion coefficients arc called splitting 
coefficients. The representation of splitting with Clebsch-Gordon coefficients rather than the commonly used 
Legendre polynomials results in an inverse problem in which each degree of differential rotation is related to a 
single splitting coefficient (eq. [56]). The combined use of the vector spherical harmonics as basis functions for 
differential rotation and the Clebsch-Gordon coefficients to represent splitting provides a unified approach to 
the forward and inverse problems of differential rotation which will greatly simplify inversion. We submit that 
the mathematical and computational simplicity of both the forward and inverse problems afforded by our 
approach argues persuasively that helioseismological investigations would be well served if the current ad hoc 
means of representing differential rotation and splitting would be replaced with the unified methods presented 
in this paper. 

Subject headings: Sun: oscillations — Sun: rotation 


1. INTRODUCTION 


An acoustic mode of oscillation of a spherically symmetric, nonrotating, adiabatic, static solar model without magnetic fields is 
typically identified by a trio of quantum numbers that represents its displacement field: w, the radial order; /, the spherical harmonic 
degree; and m, the azimuthal order of the mode. Because of the rotational symmetries of this model, the modes of oscillation are 
->l j degenerate. That is, the frequencies of the 21 + 1 modes with different m values but with the same n and / values are identical. 
These modes are said to form a multiplet. The real Sun is not so simple. Of particular relevance for this paper is the fact that the Sun 
is routing and is deforming internally so that, for example, the surface rate of rotation at the solar equator is greater than at the 

poles. This phenomenon is known as differential rotation. . 

A number of ways have been chosen to represent differential rotation mathematically. We will argue in this paper that a 
represenUtion with exceptionally nice consequences is the solar rotational velocity r to ,(r, 0, <f>), defined to be the axisymmetric 
component of general toroidal flow fields in the solar interior. A heretofore more popular, and perhaps more conceptually 
appealing, way of looking at this is that the Sun is rotating differentially and that the rotation rate fl(r, 9) is itself a function of both 
radial position and colatitude. Rotational velocity r ro ,(r, 9, 4>) and rotation rate fl(r, 9) are simply related by 

r ro ,(r, 9. 4>) = SHr, 9, <p) x r = fl(r, 0)r sin (9)d> (lt 


where Cl = Qz, z is the unit vector which points along the axis of rotation, and r is the position vector from the center of the Sun to 
position (r 9, <t>). The coordinates r = (r, 9, tj>) are spherical polar coordinates (where 6 is colatitude) and f, 6 , <j> denote unit vectors 
in the coordinate directions. However represented, solar differential rotation lifts the degeneracy of the solar acoustic or p-modes, 
splitting the frequencies of oscillation of the Sun. This phenomenon is. without doubt, the largest contributor to the splitting of solar 


Observations of split solar p-mode frequencies date from Claverie et al. (1981). The quantity and quality of new measurements 
have been steadily improving (e.g., Gough 1982; Hill, Bos. & Goode 1982; Duvall & Harvey 1984; Duvall. Harvey, & Pomerantz 
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1986: Brown 1985; Ubbrecht 1986. 1989: Brown & Morrow 1987; Tomczyk 1988; Rhodes et al. 1990). Since p-modes are 
^j'" a "l!* y H Sp it by dlfferent ' a r ° tatlon - a first-order understanding of the data requires the estimation of models of rotation that 
predict the data accurately. Furthermore, an understanding of the solar angular momentum budget requires accurate models of 
differential rotation (Gilman, Morrow. & DeLuca 1989). For these reasons, the inverse problem of infemnT^the rS Ld 
latitudinal dependence of differential rotation from observed split frequencies has been a major focus of helioseismology The reader 

rhecu^enr^r^ 0 ^^' * ‘ (1989) ’ Chnstcnsen - Dals 8 aard - Schou, & Thompson (1990), and Thompson ( 1990) for clear discussions of 
the current state of this venture. 

A necessary preface to tackling the inverse problem is the solution of the forward problem: i.e., given a model of solar differential 
rotation, determine the split p-mode frequencies. This problem was first solved by Cowling & Newing (1949) and subseauent 
treatments can be found in Ledoux (1951), Hansen. Cox. & Van Horn (1977), Gough (1981), and Brow™ SSS, 
oi the solution of the forward problem is a function of the basis functions chosen to represent differential rotation. Since differential 
rotation is axisymmetnc and even about the equatorial plane, it admits a very simple mathematical representation 

there " c *Y° Popular ways chosen to represent differential rotation, both of which are based on polynomial 
expansions of rotation rate ft rather than rotational velocity r, M . The more popular of these parameterizations is to expand Qir 8) 
in even powers of cos 8 (see Brown et al. 1 989) : expana igr, u) 


0(r, 8 ) = £ 

k * 0 , 2 , 4 ... 


fVr) cos* d 


he bu k ^ tatl ° n r f te of thc Sun and the for k > 0 describe the radial dependence of latitudinally 
Hif? CrC f nl ? r ,°, at,0n ' popularity of this representation probably derives from tradition, since it was used in early 
studies of differential rotation made from observations of the solar surface (e.g. t Howard & Harvey 1970). Thus, the use of these basis 
functions eased comparison with direct observations of differential rotation. It is generally recognized that a problem with this 
parameterization is that the basis functions are not orthogonal. This is not really an obstacle from a forward theoretic perspective 
but is troublesome mveree theoretically since future inversions for higher degree components must also redo tS lower^Lee 
?vl"rmm«^ 5lnce ' stnetly speaking, they are not independent of one another. The use of any set of orthogonal basis functions 
1988) hlS Pr ° b Cm ' In P articular ’ Legendre polynomials in cosine colatitude have seen some application (e.g^ Koraennik et al. 

G(r, 0) = f Q^/Vcos 8 ) . (3 ) 

i ,/ dtb ° Ugh the P aram ®terizations of differential rotation given by equations (2) and (3) are intuitively simple and allow straightfor- 
ard companson to other kinds of observations, their use leads to unfortunate consequences that can be entirely circumvented with 
a more judicious choice of basis functions. Problems with these basis functions include the following. (1) They do not generalize 

^Philnlv'fcR n0 H aXI n^ Ctn f flOW l {2) , They d ° n0t yicld conven,cnt ‘y to thc cl egant generalized spherical harmonic foraudira 

of Phinney & Bumdge (1973) and are, therefore, computationally cumbersome. There is a more significant problem iiiaMfi the 

way observers choose to represent splitting data as a Legendre polynomial expansion in (m//): 


+ l y m a„ P ,{m/l) . ( 4 ) 

1*1 

Thu: expansion coefficients i a„ are commonly called splitting coefficients. The use of equations (2) or (3) as basis functions for 
differentia rotation applied to splitting data represented with equation (4) leads to a serious practical problem troubling inversion 

STrUe^Iirel^fol all evenT > ° f Pr ° blCmS “ WhiCh ^ 3 * ‘ hC cstimation or Spends on Q k (r) or 

Problem 3 has been tackled in a number of wavs, including: (a) by estimating n t (r) orn,(r) for all even k < k simultaneously 
where - usu j ^* ™°™P so " l^r' est,ma ‘ ,n 8 each ^ or OiM recursively by solving first for "Sch even k’ where 
, , , f' B 7, J 989) ’ by ff™ 1 "? recombined basis functions that allow, in the high l limit, the inverse problems 

for distinct k to decouple (e.g., JCorzenmk et al. 1988); and (d) by replacing equation (4) with an alternative representation of splitting 
relative to which recombined basis functions decouple in the inverse problems for distinct k (Dumey, Hill, & Good! 

988 i y H prob !^ mS Wlth each of thcse approaches. Approaches a and b generate models that at different degrees k have 
correlated errors. In addition, approach b necessarily performs the recursion in the direction opposite from how a stableand robust 
recursive technique should be applied. A robust recursive technique would first estimate the features of the model that have the 

e r C r n ' n thc datSL In th,s case - these are ‘he longest wavelength features of the model (i.e., small k). Then these should be 
used m the estimation of shorter wavelength model components (i.e., higher k) that affect the data more subtly. Approach b does the 

^r^l th,S ^V S I ,mat, 7.L hC Sh °, rt f r A wavelength features first, it propagates errors from the more poorly constrained to thc 
better constrained feature of the model. Approach c has a limited range of applicability and approach d requires observers to 
summarize their data in a form they consider suboptimai. M 

The common problem with all of these approaches to problem 3 is that the basis functions that have been chosen to represent 
both splitting and differential rotation have been ad hoc. In this paper we take a different approach. We define the inverse problem 
explicitly m terms of the form of the solution to the forward problem. We show that there exists a natural set of basis functions with 
which to represent differential rotation such that the inverse problems for different degrees of differential rotation decouple with 
respect to data represented in the usual way (eq. [4]). Consequently, joint and recursive inversions as well as asymptotic approx- 
imations can be entirely circumvented. These basis functions are simply the vector spherical harmonic components of p rw . In 
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. pa< .. lv t0 nonaxisymmetric flow fields and yield to the formalism of Pbmney ■* 

addition, vector sphencal harmonics generau y ^ also show that the inverse problem is simplified further if splitting 

Bumdge; thus, their use also addresses problem > < d hc solution of the forward problem. These orthonor- 

«. » + . <*««*« - • ** — *« — * ■» 

represented as follows : 


= (0 -i 


I JbuPu 


(5) 


The expansion coefficients represent a set of new splitting cMfficients. havc takcn t0 determine the splitting and 

Our approach to the forward problem is motivated by the approacn ge P ^ q{ ^ Earth (e^ Dahien 

coupling of terrestrial osallations caused by ^phen^ pertu Woodhouse & Girnius 1982). An approximation that 

1968, 1969; Luh 1973, 1974; Woodhouse & DaUen 19J8 Woodho^ ivw ^ ^ ^ ^ fadial order n and harmonic 
has proven useful in terrestrial applications is t ■ th absence of asphericities, they will not be consi3ered potential 

degreej. This means that if two mote j ^ not d« iget to use degenerate perturbation theory to comp«e 
couplini partners in the presence of the asphenaty. In this ca» PP P multiplet approximation since it is accurate if the 
toe split frequencies. This approximation ,s .^ n0 ^ t ? ^^^^^^jn^frommodes aimposing other multiplets. In terrestrial 
21 + 1 modes composing a multiplet are lsolat f^ ' P ful and ^t e accura te for calculating split frequencies but is not as 
applications this approximation has proven to be hig y i accmcnts accurately, quasi-degenerate perturbation theory has 

useful for computing modal displacements. To compute r P en , f ^ only nearly degenerate with respect to the sphencal 

been used by geophysicists in which modes arc al '°7^ ' d ° J near degeneracies between modes from different multiplets satisfying 

Earth model. For the Sun. the num ber of significant^cidenta^ near degen^^ ^ ^ sQ[ar muUip i cts can accurately be 
the selection rules that govern coupling for , d te wrturbation theory can be used to compute the split frequencies. The 

X d on to ihSomXoblem ^adngof^^ 

order than first in Q/Z We will discuss this and relate them to the previously used basis 

In § 2. we defined the new basis functions forward problem for differenual rotation using 

functions given by equations (2) and (3 \ I " § 3 ’ ^ P ^ ' * trms 0 f Wigner 3 -j symbols which are straightforward to compute 
degenerate perturbation theory^The solution is expressed in terms o^ thfi wigncr 3.; symbols in terms of simple 

numerically. For low-degree differentia rota 1 , -hsch-Gordon coefficients for representing splitting data, by using equa 

polynomials in m and I. In § 4, we discuss OeveofO^ Jand are simply related to the Wigner 3 -J 

(5) as an alternative to equation (4). These c^ffi^ts fotm of thc Clcbsch -Gordon coefficients to represent spluting provides a 
symbols found in the solution of the forward P roblra ; ™ “^ n °‘ 5wcd enve the form of the inverse problems rdatmg the basis 
unified approach to the data analysis and *0 werse p ■ § ( ’ aditional sp ii t ting coefficients a,. In both formulations, a smgle 

functions for e rol both to the new splitting consents .b, udt^ can P bc re f ated t0 a linear combination of the spbtfing 

degree of the vector sphencal harroonic exp^sion h $u ^ is particularly simple, reducing to a single ‘«xn^ p o 

coefficients. However, by using the t0 the rotation rate basis functions for differential rotation. Finafly, 

6 ^ m m “ d 1 

for low-degree i coefficients are presented in the Appendix. 


( 6 ) 


. VECTOR SPHERICAL HARMONIC REPRESENTATION OF DIFFERENTIAL ROTATION 

«,*)-! ii l"lr.».«+n(r. «.♦)]. 

T h c ^ - b. f»Uy characterized h y .He .diae « vector epHencat Hanpo^c eapaneton 

coelficients ajrl »',|r). and W,(r|. ^ # ^ . m W + ^r|T, m 4) . ™ 


7>, 8. <p)= - ivyr ) r % V, Y\(8, <f >) , 


where the surface gradient operator. V„ is given by 




The function Y 


Vi=r[ V-rir-V)] = 9- + — ~ e -^- 
Y\(8, <t>) is a sphencal harmonic of degree s and azimuthal order r defined using 


(9) 


the convention of Edmonds ( 1960) as 


Y', = l — D' 


"2s + 1 U-rii T' 2 
4 it (r + silj * 


F s ( cos 6)e“ 


( 10 ) 


:'60 
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where Y~ J - i - l) r [Y, f ]\ P\ are associated Legendre functions, and the asterisk represents complex conjugation. The normal- 
ization constants in equation < 10) have been chosen such that 


I j IY[{0. &q*Y\(e, (p) sin Odd d<l) = d. t d S3 1 nn 

„o Jo 

where integration is over the unit sphere. The poioidai coefficients u\[r) and u\(r) are not independent if the anelastic condition 
[V • ipv) - 0] is imposed. 

The differential rotation velocity field r rot (r, simply corresponds to the odd-degree, zonal part of the toroidal flow field in 
equation 18) which can be written: 


where, for example: 


= - f W?(r)d a rf 4> , 

5 = 1 . 3 . 5 ,... 


r? = 


r Y° — 
1 e 1 3 


1 

2 

3 

4 



sin 0 . 



sm 0(5 cos 2 9 - I) , 


02 ) 


(13) 

(14) 




15 

16 



sin dill cos* 5—14 cos* 0 + 1) . 


05) 


The relationship between the expansion coefficients wf(r) with fl t (r) and ft,(r) can be determined by equating the representations 
in equations (I) and 1 12!. expanding each side in terms of irreducible trigonometric functions, and equating the appropropriate 
groups of terms. This procedure yields 


= 2 h r 

\ J 


fi„(r) + j n,(r) + j- 5 G 4 (r)J - 2 - 1 fl 2 (r)J 

" S(r| ■ V r [n a = 1 '’ * ih {; - ? n ‘ w ] ■ 
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. [7 n 1 


where we have truncated the sums at k = 4 in equations (2) and (3) and at s = 5 in equation (12). 


(16) 

(17) 

(18) 


3. THE FORWARD PROBLEM 

The equations governing the effect of differential rotation on solar oscillation frequencies can be presented naturally in either of 
two reference frames : a frame corotating with the average angular velocity of the Sun or an inertial frame roughly identifiable with 
the frame of observation. Since differential rotation is stationary relative to both the corotating and inertial frames, solutions to 
these equations will separate in both frames. For simplicity of use by observers, we choose the inertial frame in which to represent 
and solve the following equations. As a consequence. e ro , will be considered to include the average rotational velocity of the Sun. 

Let k -\n, [); then the displacement field for p-modes in the presence of an axisymmetnc flow field separates spatially and 
temporally : 

sir, r) = *T(r)e-;\ (19) 

We have chosen to introduce an unfortunate notational conflict and let t represent time (as well as the azimuthal order of a 
convective flow field). The radial eigenfunctions are defined as follows: 

*TM = m Ufr)Y?l8, <j>)r + 'VfrWt YT(d , <t>) , (20) 

where n L ^r) and „V^r) denote, respectively, the radial and horizontal eigenfunctions for harmonic degree / and radial order n. 
Hereafter, we drop the subscripts n and / in equation (20) and use instead V = n Uf{r\ and V = ,LJ(r). The eigenfunctions satisfy an 
orthogonality condition given by 

-s”d 3 r = rn m . m 6 n . m d n , (21) 

where p 0 is the density of the equilibrium solar model, and 

V * Po(U 2 + L}V 2 )r z dr , (22) 

Jo 

where 1} - /(/ 4- I). The scalar constant N depends on the normalization of the eigenfunctions U and V, 
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of differential rotation on p-mode frequencies. The equation of motion for the 

° ^U k ) = 0 - (23) 


„here^isa linear. self.ad.oin, spanal integro-differential operator sub,«c, » certain .deal, red boundary conditions. Tins equation 
can be rewritten using equation (191 as ,-mt 


-^f J t ) x /'o“ , i I ‘ = l) ' 

, , , ... vcinriiv field v the equation of motion is altered. In particular, the 

K2&5* iStlSSSitS “*« derivative D,'"and the dtspiaeenten, and frepuentaes of the modes 
^st be perturbed. Thus, we must make the substitutions 


(24) 


■ D t = c t + p, 


(25) 


; rot 


s k -» s k + Ss k , 
U) k — ► 0) k + d(0 , 


(26) 

(27) 


a> k -* <o k ouj , 

tvh^ <», is the degenerate frequency of the unperturbed multiple.. Lynden-Bel, dt (Writer , 19671 showed tha, the -«*» of the 

w Kn“;M “™ s ,n '* and and usiBe ' q ““°" <241 

eliminate terms, we obtain the perturbed equation of motion 


* . ^8) 
_'(<5j i ,I - lifcJkPorio. • Vs * + Ss k ~ = 0 • 

In accordance with our discussion in i 1. we will “'j™ *{£ 10 differenriLTro^ this approximation, we 

t in of L , 2,1 1 1 etgen.unc.tons of a single mult, pie, of a spherically symmetnc 

solar model as follows 


J*= I . 


(29) 


Inserting equation ,291 into equation 1281. taking the inner product of the resulting expression with jT *• integrating over the volume 
of the Sun. and using equation (241 again, we obtain 


l. and using equation again, wcuuianr _ 

j, r .. M ,^‘r - f yirf •) ' tir.tfV * j_ u.{d» 2 j P.sf • • >N‘r - Tito, j P.*T * ' ' T ’? , '‘ r } ' ° ' l3 °’ 

differential rotation : 

£ = 2to k j PoC* * • 


SSoj; 


(311 


where we have used the fact that the rtghfhand side of equn^n ,31)-“^“ ” Jr I* -*l: 

Substituting into equation (31) from equations ,20) and 1 12) and using the fact that 0o> _ 

... ^ V r . 


10? = 0) t - d(U m = to., f Z 

e = 1 .3. 


s= 1.3.5... 


. si n^sf 


(32 


where 


_ j *>/ 

;sl — 1 




(33 


i"w>l,,K ls (r)r-dr. 
Jo 


and where we have defined 


,K,tn - -p,r-!t" + 0V‘ - vw - las + W]!.* . 
HTF.-l-lT—fj l _!)• 


(34 

(35 


(3t 


F. = 


(2/ -si! 


(2/ + s + D! 


( 3 - 


The gradient operator in equation .31) acts ^ Conolri P L F +T)t “"conmbmio^^ kernel n K„. 

might loosely be called the advection , L J - L V ) and Conolis |*i, v .s,s j 
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The derivation of equation (33) is greatly simplified by use of the generalized sphencal harmonic formalism of Phinney & 
Burndge ( 19731. It is beyond the scope ot this paper to describe their formalism in any detail. However. Lavelv & Ritzwoller 1 19911 
jo describe the tormalism as it applies to the problem ot splitting caused bv a general convective flow field of which, of course 
Jifferential rotation is simply a special case. The form of equation (33). written in terms of Wigner 3-/ symbols, follows from the 
Phinney & Bumdge formalism. The product of the Wigner 3-/ symbols is simply related to the integral of three generalized spherical 
harmonics over the unit sphere (see, e.g., Edmonds 1960). It is worth noting that the appearance of the I and - 1 in the lower row of 
,he 3-/ symbol represented by H s is related to the gradient in equation (31). This reveals one of the beautiful aspects of the 
generalized spherical harmonic formalism, that gradients translate simply to index raising in the 3-j symbols. 

Another attractive feature of expressing the solution to the forward problem (eq. [33]) in terms of Wigner 3-/ symbols is the 
immediacy ot selections niies which result, in part, from properties of the 3-y symbols. Under the isolated multiple! approximation 
these selection rules are that the frequency perturbation caused by an axisymmetric flow of degree s is nonzero only if (1) the flow is 
toroidal. (2) s is odd, and (3) 0 < s < 21 As a consequence, we have written the sum in equation (32) over only odd j. Thus, under the 
isolated multiset approximation, only zonal toroidal flows with odd degree less than or equal to twice the degree of the multiplet 
contribute to the splitting. For the sake of accuracy we should point out that selection rule (1) does not derive from a property of the 
3-y symbols Rather it results from the fact that the integral kernel for a poloidal flow field under the isolated multiplet approx- 
imation is identically zero. We do not show this here since we are explicitly considering only differential rotation which is purely 
toroidal. However, this is demonstrated by Lavelv & Ritzwoller (1991). v 7 

Equation i32) together with equations (33H37) completely specify the forward problem. We have chosen to write equation (32) in 
a way that has Proven useful m geophysical applications (e.g., Ritzwoller, Masters, & Gilbert 1986, 1988; Giardini, LL and 

ood house 1988) where the coefficients ,c„ would be recognized as splitting function coefficients or as interaction coefficients. We 
argue that these coefficients. being linearly related to the model parameters are what should be estimated in anv analysis of the 

XTssrs dis ““ an a,, ' rra " , ‘ me,hod s - « * - 

ft is hoped that a major product of this paper will be formulae that are simple and efficient to use both in the forward and inverse 
problems ot differential rotation. The coefficients y* in equation (33) can be computed numerically and ail the results in this paper 
could be simply stated in terms of them. Numerical methods for computing 3-/ symbols are discussed and programs are tabulated in 
Zare (1988) However, for ease of use we will rewrite the 3 -j symbols for s < 5 in terms of polynomials in / and m by using the 
recursion relation of Schulten & Gordon < 1975). If desired, it is straightforward to extend these formulae to s > 5 by reputed 
application of the recursion relation. Setting /, = s. j 2 = , 3 = /. *, = 0. m 2 = m, and m 3 = -m in equation (5a) of Schulten & 
Gordon 1 1975) and using equations (36) and (37). we obtain 


HT.i = 


- 1 21 2s + I )mH? - s[(2/ + I) 2 - s 2 ]//?.,} . 


To initiate the recursion, the polynomial forms of the Wigner 1-j symbols appearing on the right-hand sideofequation(38)fors = 1 
arc required, and can be found in Table 2 of Edmonds ( 1960). We find by using equation (38) that for s < 5: 

H ° = l ' (39) 

»T = 2m . (40) 

T/T = 6m 2 - 11} . (4 ,) 

ff™ = 20m" — MIL} — 1 |m . (42) 

HI = 0m 4 — !0(6L 2 — 5 )m~ + 6 L}(1} — 2) , ( 43 ) 

HI = 252 m 5 - !40(2L 2 - 3)m 3 + [20L 2 (3L 2 - 10) + 48]m . (44) 

Algebraic expressions for have been tabulated in the Appendix. Using equations (39H44), y” for s < 5 can be written; 


- _ (1 
/lf~ , 

\ iff 


fjV 2 (3m 2 - L 2 M3 - Z. 2 ) 

’ 4L 2 - 3 

3 f TV 2 - 10m J j- (6L 2 - 2)m 
2 \4ir/ 4 L 2 - 3 

: 3 /_9V : ["Om 4 - KX6Z, 2 - 5)m 2 J- 6 L 2 (L 2 - 2)](L 2 - 10) 

8 \4jt/ _ (4L 2 - 3M4 U - 15) ’ 

Ji /'ll) ' : ~ 1-KX2I 2 - 3)m J + [20L 2 (3L 2 - 10) + 48]m 

16 \4ir/ (4 L 2 - 3H4L 2 - 15) 
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4 . ORTHONORMAL BASIS FUNCTIONS FOR REPRESENTING SPLITTING DATA 

Observed frequency splittings are typically represented in terms of a Legendre polynomial expansion with argument mfi % where 
the expansion coeiFicients are the splitting coefficients ,< a H as in equation (4). (For example. Libbrecht 1989 measured the splitting 
coefficients of 723 p multiplets in the range 5 < / < 60.) The sum in equation (4) has usually been truncated at M = 5 since the 
inclusion of higher orders does not significantly improve the fit to the data (Brown et al. 1989). Apparently, Legendre polynomials 
have been chosen to represent the splitting since they are well known basis functions and Legendre functions are orthogonal over 
the continuous interval [— l . 1]: 

I Pix)P,(x)dx = 6 iV . (50) 

J - 1 2i + 1 


Legendre functions are not orthonormal, but can easily be orthonormalizcd. There are two problems with their use. First, they are 
not an orthogonal basis set for representing discrete data such as split frequencies. If used to represent discrete data they are only 
approximately " orthogonal.” However, the accuracy of this approximation improves with harmonic degree / as the sampling of the 
interval [—1, 1] becomes finer. Second, and much more significantly, their use complicates the inverse problem. As we will show in 
§ 5 (eqs. [61]-[63]), by representing split frequencies with splitting coefficients n a H based on Legendre polynomials, the interaction 
coefficients „c„ are not related to a single splitting coefficient but are related to a linear combination of the „a ti . 

The purpose of this section is to point out that there exists an orthonormal set of basis functions over the discrete interval 
— / £ m <, /, and that these functions are simply related to the yJJ functions. These basis functions are Clebsch-Gordon coefficients. 
Furthermore, in § 5, we show that the use of Clebsch-Gordon coefficients to represent splitting further simplifies the inverse problem 
(eq. [56]). 

Inspection of equation (32) suggests that a natural, alternative way to represent splitting caused by differential rotation would be 
to use the functions vJJ as basis functions rather than Legendre polynomials. For fixed /, the yJJ functions are orthogonal in s. Their 
orthogonality relation can be deduced from the orthogonality property of Wigner 3 -j symbols given by equation (3.7.8) of Edmonds 
(I960). We rewrite this for our purposes as 



Combining this with equations (33) and (36), we obtain the orthogonality relation for the y" functions: 


where we have defined 



y?i y?i 


g, g, 

2s + 1 F, F r 1 


G, * y?JH ? . 


(52) 

(53) 


We note that G, is independent of m since the m-depcndencc of y" is given by as can be seen in equation (33). 

Since the terms on the right-hand side of equation (52) can vary widely in size with s, especially for high /, theyJJ should not be 
used as basis functions for the splitting However, equations (51) and (52) suggest a set of basis functions which are orthonormal on 
the discrete interval - / < m <, /. These orthonormal functions are simply related to the yJJ functions as follows: 


where 


am _ 


= (2s + I) l ' 2 F 5 tf7 


I nfift = 4, . 


(54) 


(55) 


Equation (54) can be evaluated explicitly in terms of polynomials in m and / for s < 5 by use of equations (39H^ anc * (37). Higher 
degree can be computed by using the formulae provided in the Appendix. By equation (3.7.3) of Edmonds (1960), it can be seen 
that the fi" coefficients are simply Clebsch-Gordon coefficients. Equation (5) would then be used to represent splitting. 


5. THE INVERSE PROBLEM 

Once the new splitting coefficients n b u have been estimated using equation (5), the inverse problem for differential rotation can be 
reconstructed immediately since the n b 3l coefficients are simply related to the n c si coefficients as follows: 





uf(r)/C s (r)r 2 dr 


(47T) W2 

(2/4- I) L 2 H l ,F, " 


(56) 


where we have employed the notational simplification c s = n c zi and 6, = n b zl . Given estimates of the interaction coefficients c„ 
equation (56) defines a linear inverse problem for differential rotation. 


;64 
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!th0l |, Sh CStimating j hC ? eWSpil , tting coetTicients b , w,th the orthonormal basis functions 0" is more stable than current methods 
and results in an exceedingly simple inverse problem leq. [56]). some observers mav choose to retain the wav splitting data have 
traditionally been represented. Thereiore. we also seek expressions which relate the estimated splitting coefficients a to the 
interaction coefficients ,c„. W, do I his only for , . 1. y 5 here, although the method we entpiov can S uL“ toTSiS tagto 
degree results. First, identify equations (32) and (4): ' prouucc mgncr 

5 

1 £ a,P,(m/[) = V G, H?c , , 


S= 1 . 3.5 


157) 


where we have set a, = and we have defined G, in equation (53). Then, simply equating terms with the same odd power of « 


a > " : a 3 + ¥<3, = 2 C,c, + (4 - 12 L 2 )G 3 c 3 + [20L 2 (3L 2 - 10) + 48]G,c 5 , 
a 3 - ia, = 8/ 2 G 3 c 3 - 56/ 2 (2L 2 - 3)G 5 c, , 
a 5 = 32 l*G,c s . 


(58) 

(59) 

(60) 


for 75 c45H4,]) “ d * cw] - [44]) in, ° ,h ' •*«■■■ »» b. 


w*{r)K Ar)r~ dr = 2 /- 
V 3 


.0 

■* R - 


-3 = I n?(r)JC 3 <r)r^</r = - /- 

Jo V 7 

»**© 


~ I a 3 (, l\ a 5 /\ 7 27 9 

. ' ’ r 2i V /,/ + zi ( 3 + :/ ~ 4 1 1 ^ 4/ j 

+ 3) r a, I 21 y 

— - 37 ( 7 - 37jj 


?r (2/ — 1X2/ 


= 1 w?(r)K j (r,r 2 dr= /— — ~ 3)(2/ ~ 1 H + 3 X2/ ± 5) 

Jo ‘ V 11 15/ 4 Lfl5j - 


(61) 

(62) 

(63) 


q a ions (61 H63) constitute three independent inverse problems for the radial functions w?(r), vv^fr), and w°(r). Equation (63) 
follows immediately from equation (60) Equation (62) was obtained by substituting equation (63) into equation (59) and solving for 
1 ^ Equation (61) was obtained by substituting equations (63) and (62) into equation (58) and solving for c,. For clarity, it should be 
pointed out that the number of terms on the right-hand sides of equations (61H63) depends on the accuracy with which the split 
frequencies are estimated. As frequency estimates become more accurate, the number of terms will increase. 

r T™ S SCparat J on ,nt0 threc independent in verse problems, each uniquely identified with a single harmonic degree j of differential 
tation, has been made possible by use of the vector spherical harmonic basis functions given by equations (6) and (12). The 
corresponding inverse problems for Cl, and n k do not separate so nicely. The advantage of using the Clebsch-Gordon coefficients as 
basis functions for splitting is readily apparent by comparing equation (56) with equations (61H63). The use of the Clebsch-Gordon 
coefficients reduces the linear combination on the right-hand side of equations (61H63) to a single term in equation (56). 
rJn 15 be J° nd , '7? "FI ?/ thlS papCr ‘° d ' sc * 2ssthe lar 8 c and well-known variety of approaches to linear inverse problems. The 
lomi p l rCC ,m77i ^ 8 P A P ? rS Wh ' ch d,scuss a PP roachcs t0 these problems in some detail: Backus & Gilbert (1967, 1968, 

JJJ , ( J™ten*en-Dalsgaard et al. (1990). Once *v?(r), w°(r), and w?(r) have been estimated from equation (56) or 

from equations (61H«) by whatever inverse method has been chosen. 6, and Cl k for k = (0, 2, 4) can be computed a posteriori, if 
desired, by equations U6H 18); 


Q 0 (r) = — ==- [ v 3w°(r) + v 7w°(r) + v llw°(r)] , 
v 4irr 

Q :( r ) = “==- [ v M(r) + v TTvvj(r)] , 

v 4nr 

n*(r) = “ 4 - [v Hwf(r)] , 

^ 0 (r) = [ v '3 w°(r) - 4 V 7w°(r) + ¥ v 'TTw?(r)] , 

\ 

[ v 7w?(r) - i v T7w°] , 

2 v 4nr 

= 7^ Cv TTw?(r)] . 

8 ^/ 4tt r 


(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 
(69) 
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From equations (61H63I. it can be easily seen that in the limit of large 1. the expansion coefficients w?<n. w’ftrt. and w?(r) 
oniv on a,. u 3 , and a., respectively. Thus, in the large i limit, these equations can be rewritten: 


\ l 2 


dr . 


€ 


-tUtt) 


d « . 



6. SUMMARY AND CONCLUSIONS 

At its incemion the motivation for this paper was that there did not appear to have emerged a fully coherent, unified treatment of 
the^fo^Tr^Tnd^nverse^roblems for differential rotation. As a result a number of practical problems have beset mve^om for 
inferential rotation In this paper, we present a general formulation of both the forward and inverse problems for differential 
“ welf aT specific formulae useful to observers with the intent of simplifying solar differential rotation inversion. There are 
two main Doints of the paper ( 1) If differential rotation is represented by rotational velocity, defined as the zonal, odd degree parto 
rS-Ef ESSSi decomposition of a general convective field in the solar interior, rather than the conunotfiyus^ 
rotation rate represented with ad hoc basis functions, then several significant problems currently facing inversions for differential 
rotation disappear. In particular, the inverse problems for different degrees of rotational velocity are linear and decouple, so that 
^dependent inversions for each decree of structure can be performed without approximation. (2) The inverse problems a 
'icmficantlv simplified further if Clebsch-Gordon coefficients are used as basis functions to represent splitting. Thc ^^ b '£i r ®" 
basis functions are genuinely orthonormal on the discrete interval -l <m< l and regression matrices comprising 
optimally well-conditioned. As a consequence, we highly recommend that the vector spherical harmonic representation of ro |* * 
and the Clebsch-Gordon coefficient representation of splitting be adopted to replace the ad htK rep^enUtion^ 
employed^heretofore. We have presented formulae relating the interaction coefficients c, to both the new splitting coefficients b, as 

^^eophyfica^^'^ne^cThas^hown that it ’is useful to tabulate splitting data in terms of the splitting function or inaction 
coefficients c sincrihese coefficients are linearly related to the structures producing the splitting. However, in the Sun, there 
number of different kinds of axisvmmetric asphencities that could produce splitting. In addition to differential rotatiom p OTomena 
which have been discussed as potentially large enough to affect p-mode frequencies measurably include lateral density varwuons 

variations, asphenat.es in the figure of the Sun, and large-scale, dominandy quadru^lar 
magncu^fields Each of these mechanisms affects splitting differently. However, with the possible exception °* **^ ?^ 
component of magnetic fields (D. Gough, 1990, personal communication), in each case ^ 

guarantees that the solution to the forward problem can be wntten in a form identical to that of equation (3-), but with 

dHclmVin Si for mch mgchinism of spiriting, Thus. an inspection oi equation ,33) reveal, that ih. Mtouton 
m .. , shoul d no i be tabulated since thev difTer in detail among the various sources of splitting However, rite r3 coefficients 
mme ,a™;r»m«“ ofTpl ri,ng are similar ,» that each is linearly related to the same Clebsch-Gordon eoeffiemot Consequently. 
« !SSt new splitting coefficients », estimated relative to the Oelrich-Gordo, coeffiomtt tepm rimtaj.oo of 
-nlittinc For each of the sources of splitting mentioned here, these coefficients are linearly and simply related to th pp p 
Ssls fuLfi^s differential rotation we have shown that the appropriate basis fane* onsare | pven 

by eq. [12].) Of course, if multiple sources are causing splitting, then models of the various mechanisms woul avc o 

^h™uchoin the naoer. we have attempted to present formulae that would prove to be easy to use. In conclusion, it « worth 
resenting a brief review of the most important of these. The differential rotation basis function p ro , is defined in equations 
H2) The general solution to the forward problem is given by equation (32) and the equations tmmed.ately 
form of the forward problem for s < 5 is presented in equations (45H49) where polynomials in m. I, and s replace & * 

symbols of the general solution. [Formulae useful for calculating the solution to the forward problem analytically forb J^ d g L 
symbol * ™ he Equatlon (3 4) is the basis for the inverse problem with respect to the «»tots.The 

cfrthonormai a basis°funcdo^ I Clebsch-Gordon coefficients, to represent splitting are defined in equation (54), ^ tabffiat^ .n ffie 
Annendbtandffie inverse problem with these coefficients is given by equation (56). Relative to the coeffiaen^ the = 
probes are pven S equations (61H63). Once r rol has been estimated, the rotation rate can be computed. ,f desired, with equations 

(64H69). 
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We have argued in this paper for the following representation of splitting data: 

M 

'■jZ = CZJ„ t V „b,JT, , 

5 = I 

■' here are the new splitting coeflicients and are Clebsch-Gordon coeflicients defined as 


(Al) 


K = <-5 + d 12 f,h: , 


(2/ -si! T ;2 
(2/ + s + l)!_ 


(A2) 

(A3) 


E A X a P iL e aiiw f0r ^ r l 0r , 1 5 uT ^ CD by eqUaUons (40H44) in thc bod y of the tcxt and for 6 < s < 11 below by equations 
A4M A9). We note that with the H, coefficients provided below, it is straightforward to calculate the yJJ coefficients for a > 5 by use 
or equation (JJ). ' 

The //" coeflicients for 6 < s < 1 1 are given by 


tf" = 924m 5 - 420m 4 (3L 2 - 7) + 84m 2 ! 5 L 4 - 25 L 1 + 14) - 20 L 2 (L 4 - 8 1} + 12) , 

H - = 3432m 7 - !848m s (3L 2 - 10) + 168m J (15L 4 - 105 L 2 + 101) - 8m(35L 5 - 385L 4 + 882L J - 180) , 
H” = 12870m s - 12012m 6 |2L J - 4) + 2310m 4 (6L 4 - 56L 2 + 81) 

— !2m-|210L° - 3045L 4 - 9898Z-* - 4566) + 70L 2 (L° - 20L* - 108E 2 - 144) , 

= 48620m 9 - 17160m 16 L 2 - 35) * l2012m 5 (6L 4 - Ills - 145) - 440m 2 

< i42E® - 777L 4 + 3402L 2 - 2630) + 12mtl05L B - 2660Z? + 18844L 4 - 36528L 2 + 6720) , 

H "o = 184756m 10 - 145860m 8 (3E 2 - 22) + 12012m 5 <30L 4 - 450L 2 + 1199) - 2860m 4 

x|42L 6 - 966 L 4 + 5481L 2 - 6248) + l32m 2 (105L 8 - 3290L 5 + 29680L 4 - "8900L 2 + 32208) 

— 252L 2 (L 8 - 40L 5 + 508L 4 - 2304L 2 + 2880) , 

H^x = 705432m" - 1847560m 9 (L 2 - 9) + 58344m'(30L 2 — 550 1} + 1869) — 120120m 5 

x| 6 L 5 — 168L 4 + 1 199Z. 2 — 1873) + 1 144m J (105L 8 — 3990i? + 44730Z? — 156200L 2 + 105228) 

— 24mf 23 1 L 1 0 - 11165L 8 + 174328L 5 - 1006764L 4 + I771440L 2 - 302400) . 


(A4) 

(A5) 

(A6) 

(A7) 

(A8) 

(A9) 
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ABSTRACT 

In r his paper we derive a theory, based on quasi-degenerate perturbation theory, that governs the effect of global- 
scale. steady-state convection and associated static asphericities in the elastic-gravitational variables (adiabatic bulk 
modulus ac. density p. and gravitational potential o) on helioseismic eigeni'requencies and eigenfunctions and present 
t formalism with which this theory can be applied computationally. Eigenfunctions are simply seismic displacement 
patterns. The theory rests on three formal assumptions: (Ij that convection is temporally steady in a frame corotating 
with the Sun. i2) that accurate eigeni'requencies and eigenfunctions can be determined by retaining terms in the 
seismically perturbed equations of motion only to first order in p-mode displacement, and (3) that we are justified in 
retaining terms only to first-order in convective velocity as well (tiiis assumption is tantamount to the requirement 
that the convective flow is anelastic). The most physically unrealistic assumption is f 1), and we view die results of this 
paper as the first steps toward a more general theory governing the seismic effects of time-varying fields. Although 
the theory does not govern the seismic effects of nonstationarv flows, it can be used to approximate the effects of 
unsteady flows on the seismic wavefield if the flow is varying smoothly in time. The theory does not attempt to 
model seismic modal amplitudes since these are governed, in part, by the exchange of energy between convection and 
acoustic motions which is not a part of this theory. 

The basic reference model that will be perturbed by rotation, convection, structural asphericities, and seismic 
oscillations is a spherically symmetric, non-rotating, non-magnetic. isotropic, static solar model that, when subject to 
seismic oscillations, oscillates adiabaticaily. We call this the SXRXMAIS model. An acoustic mode of the SXRNMAIS 
model is denoted by k = ( n . /, ml, where n is the radial order. / is the harmonic degree, and m is the azimuthal order 
of the mode. 

The main result of the paper is the general matrix element for steady-state convection satisfying the 

anelastic condition with static structural asphericities. It is written in terms of the radial, scalar eigenfunctions of 
the SXRNMAIS model, resulting in equations (92)-(112). We prove Rayleigh's principle in our derivation of quasi- 
degenerate perturbation theory, which, as a by-product yields the the general matrix element. Within this perturbative 
method, modes need not be exactly degenerate in the SNRXMAIS model in order to couple, oniy nearly so. General 
matrix elements compose the Hermilian supermatrix Zk'k ■ The eigenvalues of the supermatrix are the eigenfrequency 
perturbations of the convecting, aspherical model and the eigenvector components of Z^k are the expansion coefficients 
in the linear combination forming the perturbed eigenfunctions in which the eigenfunctions of the SXRXMAIS model 
act as basis functions. We show how helioseismic synthetic seismograms can be computed using the supermatrix. 

The properties of the Wigner 3-j symbols and the reduced matrix elements composing produce selection 

rules governing the coupling of SXRXMAIS modes that hold even for time- varying flows. We state selection rules for 
both quasi-degenerate and degenerate perturbations theories. For example, within degenerate perturbation theory, 
only odd-degree s toroidal flows and even degree structural asphericities, both with s < 21, will affect the splitting and 
coupling of acoustic modes with harmonic degree l. In addition, the frequency perturbations caused by a toroidal flow 
display odd symmetry with respect to the degenerate frequency when plotted from the minimum to the maximum 
perturbation. ^ 

We consider the special case of differential rotation, the odd-degree, axisymmetric, toroidal component of general 
convection, and present the general matrix element and selection rules under quasi-degenerate perturbation theory. 
We argue that due to the spacing of modes that satisfy the selection rules, quasi-degenerate coupling can, for all 
practical purposes, be neglected in modeling the effect of low-degree differential rotation on helioseismic data. In 
effect, modes that can couple through differential rotation are too far separated in frequency to couple strongly. This 
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,5 not the case when we consider nonaxisymmetnc flows ana asphericitics. In ihis case, near degeneracies will reiruiariv 
’.'cur. and couplings can be Hativeiv <uong -specially amons SNRNMAIS modes within the same multiplet. 

All derivations are performed ana ail -Gluttons arc presented in the ccrotating frame. Equation { 1$) shows how to 
* ran storm the eigenfrequencies and eigenfunctions iti.tiie corocating frame into an inertial frame. The transformation 
has the effect that each eigenfunction in the inertial frame is itself time-varying. That is. a mode of oscillation, which 
.s Termed to have a single frequency in the corotating frame, becomes multiply periodic in the inertial frame. 

I. INTRODUCTION 

Helioseismic images of the acoustic velocity held of the Sun are providing new and unique information about solar 
structure and dynamics. To continue to utilize effectively the information provided by the continually improving data 
fees wul require a thorough understanding of the way in which a number of solar structures and processes affect 
neiioseismic data. It is upon such an understanding of these forward problems that any future inversions will rest. 

We consider here the helioseismic effect of one such solar process: convection. In particular, the purpose of this 
paper is to present a theory that governs the effect of large-scale, steady-state convection, with associated aspheri cities 
:n the structural elastic-gravitational variables (adiabatic bulk modulus k. density p, and gravitational potential o ) t 
on helioseismic oscillations. Many studies have been completed concerning differential rotation, the long- wavelength 
ixisymmetric component of convection I e.g., Duvall and Harvey 1084: Brown 1985: Duvall et al. 1986: Libbrecht 
A'>6. 1989: Brown and Morrow 1987: Rhodes zt al. 1990: Brown t it al. 1989: Thompson 1990: and Ritzwoiler and 
Laveiy 1991V However, to date, studies of the seismic effects of non-axisymmetnc convection are rather sparse. In an 
.symptotic treatment. Gough ana Toomre [ 1983} calculated the frequency shift of an acoustic mode due to advection 
by a purely horizontal how. Brown i 1984) calculated the influence of turbulent convection on modal degenerate 
frequencies. The scattering of sound by an isolated, steady laminar compact vortex was considered by Bogdan (1989). 
Hill (1983. 1988, 1989) has used a ray-theoreric method and has attempted to infer horizontal convective velocities 
near the solar surface using helioseismic data. All of these studies make restrictive assumptions about the geometry 
of the flow field including either that the flows are horizontal in a plane-parallel medium or demonstrate cylindrical 
symmetry, and none attempts to model wavefront distortion and deflection caused by convection. In summary, to the 
best of our knowledge no general theory for rhe effect of convection on helioseismic oscillations currently exists. 

The theory presented in this paper differs from these previous studies in the following ways. (1) Our theory is 
nonasymptotic. In principle, the results are accurate for ail wavelengths and frequencies of helioseismic oscillation. (2) 
It is derived within a spherical geometry. Previous investigations that modeled convective effects within a nonspherical 
geometry are appropriate for short- wavelength convection but inappropriate for the largest scales of convection which 
are the subject of this paper. (3) The theory presented here makes no assumptions about the geometry of the flow. We 
represent general non-axisymmetric flow fields comprising both poloidal and toroidal components in terms of vector 
spherical harmonics, which are complete basis functions for a vector field in a sphere. (4) Our approach is modal- 
theoretic rather than rav-theoretic. From a traveling wave perspective this means that both wavefront deformation 
as well as the perturbation in local sound speed by convection are modelled. In modal-theoretic language, convection 
results in modal coupling as well as splitting. 

In a later paper we will implement the theory presented in this paper using a numerical simulation of large-scale 
convection and discuss the observational consequences of the theory. In particular, we will show that the helioseismic 
frequencies, displacement patterns, and line-widths of an aspherical solar model are appreciably altered relative to 
the corresponding quantities calculated from a model with differential rotation alone. 

a I Modal Xotation and Terminology 

The basic reference model to which all subsequent structural perturbations and processes will be added is a solar 
model that is spherically symmetric, nonrotating, nonmagnetic, isotropic, and static, subject to adiabatic acoustic 
oscillations. We refer to this as the SNRNMAIS solar model. An acoustic mode of oscillation of any solar model is 
defined to be a characteristic spatial displacement pattern that oscillates with a single frequency. 

An acoustic mode of a SNRNMAIS model is uniquely identified by a single triple of quantum numbers (n,l,m) 
that denote, respectively, the radial order, harmonic degree, and azimuthal order of the mode. A modal frequency 
for such a model is simply the degenerate frequency of the multiplet n 5f that comprises the (2/ + 1) modes with 
identical n and / values. Any symmetry-breaking agent such as rotation, magnetic fields, or convection will lift this 
(2/ -r l) degeneracy and split the frequencies of the modes composing the multiplet. We call any model with such a 
symmetric-breaking perturbation a non-SNRNMAIS model. A major goal of this paper is to provide formulae with 
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v/nich co calculate the modal eigenfunctions and eigenfrequencies of both the SNRNMAIS and non-SNRXMAIS solar 
moueis. As we snail later see. t he perturbations of the non-SNRXMAIS ^oiar model will be assumed to have small 
masnmuae and be stationary 111 a frame coroutine with the Sun. If the symmetry-breaking agent is axisynimetri<\ 
-.s t> diilerenuai rotation, then to a good approximation the spatiai structure of each mode will remain specified by 
*:ie same triplet of quantum numbers. For a general, tion-axisymmetric perturbation such as the convective fields \\>* 
;:nsider here, the perturbed eigenfunction »or spatial displacement pattern) of each mode is a linear combination of 
'lie eieenfunctions of the 5NRNMA1S moael. We call this phenomenon oscillation-oscillation coupling or interaction 
-o distinguish it from oscillation-convection coupling, the exchange of energy between seismic and convective modes. 
The acoustic modes chat are said to couple as a resuit of a convective flow or an asphericitv in the elastic-gravitational 
vanades are 5NRNMAIS modes. The modes of the non-SNRNMAIS solar model do not couple. 

Two modes of the SNRNMAIS soiar model whose spatial eigenfunctions are orthogonal are said to be isolated 
from one another. Two or more modes tnat are not isolated from one another can couple when the reference model is 
perturbed either by a structural perturbation or a convective flow. A multiplet composed of modes whose combined 
eizenspace is orthogonal to the combined eigenspace of the modes composing all other multiplets is said to be isolated 
or seif-coupled. The degree of coupling between SNRNMAIS modes is a function of a number of factors, among which 
--re the strength of the asphericitv or convective field producing the coupling, the proximity of the eigenfrequencies 
of the modes, the relation between the geometries of the perturbing field and the oscillations which is encoded in a 
<rt of analytical angular selection ruies. and the similarity of the radial eigenfunctions of the two modes. When two 
SNRNMAIS modes k - in. /.mi and h r = iri.i'.m') couple, the strength of interaction is described by the general 
matrix element . . The matrix n n composed of ail the general matrix elements ror the multiplets n Si and 

is of dimension 1.2/' — 1) x \2! — I) and is called the general matrix. The square general matrix H nn n is called the 
splitting mat tlx and governs self-coupling. The eigenfrequencies of non-isolated modes that couple within or across n 
or / are the eigenvalues of an assemblage of block diagonal splitting matrices and off-block diagonal general matrices. 
The enure assemblage is called the supermatrix Z*';.. 

Since the seismic modes of the SNRNMAIS solar model are spheroidal, their spatial vector eigenfunctions st(r) 
-or displacement patterns) may be wmten in the form 






(I) 


where n f7j(r) and n^/(r) are the scalar radial eigenfunctions for harmonic degree l and radial order n . With the 
gravitational potential scalar eigenfunction. a O|(r), and its radial derivative, a th( r ) and n^/(r) and their radial 
derivatives form the set of scalar radial eigenfunctions. The coordinates (r.0,d) are spherical polar coordinates 
(where 9 is coiatitude) and r, 9.o denote unit vectors in the coordinate directions. The surface gradient operator is 
given bv 

^ j = r(V - r(r ■ V)). (2) 


The function is a spherical harmonic of degree l and azimuthal order m defined using the convention of Edmonds 

i 1060): 


r 1 

I W7 {9, o)] m Y™{(h o)s\n9d9do = Sm'mtn 
J a 


(3) 


where integration is over the unit sphere. Henceforth, we drop the subscripts n and l in equation (1) and use instead 
U i-i ( r ) , U f tV(r) and so on. The SNRNMAIS spatial vector eigenfunctions satisfy an orthogonality condition 
£iven by 

H) 


where 


j ‘ T — N ^rn'm^n'n &i f l 

-V = [** p Q [UU‘ +t(l+ l)VV’]r 2 dr. 
Jo 


( 3 ) 


and d 3 r = r : sin OdOdodr. Henceforth, an integral sign without limits, as in equation (4), will denote a three- 
dimensional integration over the volume of the solar model. The scalar normalization constant .V depends on the 
normalization convention of the eigenfunctions U and V. 

Perturbation theoretic techniques are usually employed to calculate split acoustic mode eigenfrequencies and 
perturbed eigenfunctions. In this paper, we will show how quasi-degenerate perturbation theory can be applied to 
determine these quantities for a convecting solar model with associated aspherical perturbations in the structural 
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-^asuc-sravnationai parameters * iesisiry t > and adiabatic buik modulus k. In particular, we will use the eigenfunctions 
:fr he FNRNMAIS model as b:ws functions ro represent the perturbed eigenfunctions sar .t): 


. L it a: J 

W~ will snow now ro determine the appropriate eigenspace /V required to represent the perturbed eigenfunction, will 
ierive rhe expansion coefficient a\ for ^acn component of the eigenspace to calculate the perturbed eigenfunction, and 
wui derive expressions for the spiit eigenfrequency of the mode = w rf / t dw : *. 

The major theoretical resuit of the paper is analytical expressions for the general matrix elements that compose 
the supermatrix (or splitting matrix in the rase of seif-coupling). The perturbed modal frequency is simply 
m eigenvalue of the supermairix (or splitting matrix), and the expansion coefficients are simply the eigenvector 
components a\. We wish to emphasize at this point that we will not attempt to present a theory that accurately 
predicts modal amplitudes, but only modal eigen frequencies and perturbed eigenfunctions (or displacement patterns). 
The formal assumptions of the theory discussed in §I.b will reflect this point. 


b) Assumptions and their Implications 


Although the theory we present in this paper is more nenerai than previous work, its application is restricted both by 
Tactical considerations and by riie set of assumptions upon which it is formally based. The major practical limitation 
.2 mat the convective structures considered mould be global in extent. For example, although it is theoretically 
possible to represent a single smaii-scaie convective vortex in terms of vector spherical harmonics, there are better 
representations and doing so would probably be a misuse of this theory. Thus, though the theory holds lor all but very 
short wavelength, turbulent convection, it will be most usefully applied to long wavelength flows. There is a caveat: 
spatially repetitive small-scale structures, such as the solar granulation, can be well represented by vector spherical 
harmonics and are not beyond the practical limitations of this theory. Another practical limitation is that the theory 
:s not particularly useful for ^-modes since they have very small amplitudes in the convection zone. Therefore, though 
'he theory governs the effect of convection on //-modes our discussion will center on p-modes. 

Much more restrictive are the following set of formal assumptions. (1) The convection is steady in time. As we 
will discuss in §11, this assumption is necessary for the perturbed equations of motion to separate. The asphencities 
in the structural elastic-gravitational variables will also be assumed to be time invariant. (2) We retain terms in the 
seismically perturbed equations of motion only to first-order in p-mode displacement and quantities that depend on 
it. Thus, we derive and use linearized equations of motion. (3) We also retain terms in the seismically perturbed 
equations of motion only to first-order in convective velocity. This is done so that acoustic oscillations and convection 
do not exchange energy and to this extent can be considered independently. This is tantamount to the requirement 
that the convective flow field is anelastic. We wiil discuss briefly the implications of each of these assumptions in turn. 
Arguments are presented to justify assumptions (2) and (3) in §I.e. 

i 1) If convection is steady in time, each identically directed acoustic wave that propagates through a given region 
will experience the same convective elTect. In particular, multiply orbiting waves propagating along near great-circles 
will experience a constructively accumulating erFect in that region. In this case, the split modal frequency associated 
with the propagating wave wiil be rime invariant. If the convective state changes appreciably during the time it takes 
an acoustic wave to execute a single orbit, then the convective effect will vary between orbits. In fact, the effect may 
destructively accumulate. Consequently, modal frequencies would be time varying, leading to an effective spectral 
iine-broadenins. This line-broadening is not a part of the theory presented in this paper and the seismic effect of 
aspects of convection that are rapidly evolving in time cannot be determined from the results presented here. Of 
particular significance is the fact that the elfect of the shearing of sectoral or banana cell modes of convection by 
differential rotation cannot be modeled within this theory. Rather, the results in this paper represent the first steps 
toward constructing a more general theory that governs time-varying fields. 

Although the results in tins paper are correctly applied only to steady-state convection, they may be most use- 
ful if seen to provide instantaneous frequencies and displacement patterns for a time-varying convective field. These 
instantaneous frequencies would be accurate over the lifetime of the convection cell which, for long-lived inodes of con- 
vection. may he appreciable. In this case, the steady-state assumption would amount to a short-time approximation. 
For example, since the shearing of convective patterns takes time to develop, the results presented here are applicable 
until the shearing effects accumulate. The numerical simulations of Glatzmaier and Gilman ( 1981. 1982) show that 
some components of flow have lifetimes on the order of weeks. Furthermore, there are certain observable solar features. 
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;:i particular active loii£itudes ana coronal hoies. that appear to evoive relatively unsneared by differential rotation, 
rf rhese features are somehow anchored at depth in convective structures. - 1 lion their existence is turther evidence tur 
a relatively stable component of ilow ueep m the convection /one. Evidence tor the existence ot solar Aiant ceils is 

■iscussed m U.d.j . 

From the view of the ste,iwy--tate assumption as a short-time approximation, it is straightforward to implement a 

numerical formalism to approximate the rime-varying seismic waveheld if we asume that the variations in convection 
ire temporally smooch. We would calculate a tune sequence of instantaneously valid eigenirequencies and eigenvectors 
-n a coarse set of time knots where at .-ach knot the how held is assumed to be stationary. We would then interpolate 
rue eisenfrequenctes ana eigenvectors onto a liner time grid and allow the waveheld to evolve coniinuouslv through 
-ach of the intervals between the knots. 

1 * 2 ) Xesiecting higher order terms than first in the seismically perturbed quantities amounts to neglecting seismic 
seif-advective eifects evidenced through acoustic three-mode coupling. In particular, the self-advection of the dis- 
placement field is neglected which is tantamount to assuming that the total seismic displacement in a region is much 
smaller than the displacements produced by convection during the passage of a wave. (We also neglect all source 
-erms such as entropy and internal energy fluctuations caused by the oscillations.) As we will discuss in §§I.c an .e, 
:he accuracy of this assumption improves with depth. The application ot the theory will be most accurate for seismic 
paths below the strongly -uper-adiabatic layer near the solar surface where turbulence is most \igorous. 

3 ) \\> nedect ail terms second-order in the convective velocity. There are two mam types of second-order terms 
mat we discard, auvective terms ana Reynold's stresses. Discarding the former amounts to assuming that convective 
roomies are relatively -mail. Ignoring Reynold's stresses, which are proportional to the Lapiacian ot the con\ecti\e 
vejocitv. requires that convective wavelengths be relatively large, and implies that convection-oscillation coup mg 
.5 nesiiected so that there is no mechanism by which convection and the acoustic oscillations can excnange energy. 
In particular, we assume that convective rlows do not generate seismic waves and. therefore, we require that the 
flows satisfy the unperturbed continuity equation commonly called the aneiastic condition. This condition eliminate* 
potential sources, sinks, and cavitation in the flow field. Thus, we view convection as a sort of passive background oil 
which acoustic oscillations are superposed. It deforms seismic wavefronts and perturbs local sound speeds, but does 
not exchange energy with acoustic waves. The assumption that second-order terms tn convective velocity and t le 
Reynold's stresses can be ignored is poor near the surface but. as with formal assumption (2). improves wit ept i 


below the photosphere. . , , .. ,, 

In summary, the implications of these assumptions are that the convective fields to which t e t eon is app ica e 
should be global in extent, relatively iong wavelength, steady in time or at least relatively long-lived, and well below 
the photosphere. Giant-cell convection satisfies these criteria and provides the best target for the application of the 
theorv presented herein. In the remainder of this section we will discuss solar convection, review the evidence for the 
existence of giant-cell convection, and attempt to justify the use of linearized equations of motion to determine the 
seismic effect of giant-cells. 


cj Solar Convection and jts Seismic Effects 


Observation of the distinct cellular motions of granules and supergranuies suggests that there are preferred sca es 
of motion for thermal convection. The common picture of convection is that the Sun contains a muitip icity 0 sca es 
of motion ranging from the Kolmogoroff microscales at the short end to differential rotation which is global in extent. 
At intermediate length scales, convective modes are thought to be organized into granules, supergranules, giant cel Is, 
and energy-bearing eddies. Temporal scales also range from a few minutes for granule overturn times to weeks for 
the largest scale of giant cells deep in the convection zone. Goldreich and Kumar (1988) present a recent review 
of turbulence. Bray et at. (1984) and Gilman (1987) provide overviews of the physics and morphology of granules. 

supergranules, and giant cells. . , 

In order to discuss qualitatively the likely general characteristics of convection below the photosphere, we loo • 
to mixing-length theory for guidance. In the mixing-length picture of convection one would take the mixing length, 
Macii number, and the time and velocity scales of convection to be given, respectively, by 

H - q H r 


U ~ 


tjF,HQp 


1/2 


4 kV- c? T 


1/3 


v H ~ c.U 


( 3 ) 

( 0 ) 


Of POOS QUAU'-Y 


vnere = ? ! {pq) :s - ;ie pr^sure scale h-'u-it. * is the ratio of : he mixing length to f/„. p is pressure, F is 
■ mp-nrure. j is tne gravity. i\u: >pec:ric beat at constant pressure, and c ;s tne souna speed. We have set 
= ; 1 - where ;s the ratio of the nos pressure to the total pressure. The convective riux can be calculated 



vnere £ is the solar luminosity. W* have u-ed equation 1 14.53) of Cox and Giuii (1963) ro obtain equation (3). 
Figure i is a plot of the 'haractensnc length, velocity, and time scales of convection predicted by equations (7). (9), 
and *10) using: the solar model of Podsiaulowski t 1939) with a taken to be 1.305. The predicted time and velocity 
-caies near the surface correspond well with observations of solar granulation. 

Convection at all depths m the convention zone will affect helioseismic oscillations. The soiar p-modes have scales 
'iiai range in size from tile smallest to the largest of the convective motions and the dominant modal frequencies 
'cinciae with the characteristic overturn times for convective motions near the soiar surface. Since the energy and 
‘he characteristic iength and time scales of convection vary with depth, the physics of interaction between acoustic 
modes and convection wiil necessarily also vary. For example, granule and sub-granule scale motions are thought to 
r.e me source of the acoustic osciiiations iGoidreich ana Kumar 1933). To model the total effect of convection at this 
-•*aie on tne acoustic oscular ions wound be very difficult as it would involve modelinz convection-oscillation coupling 
:: audition to osciiiation-jsculaticn coupling. As the formal assumptions indicate, we have set for ourselves a simpler 
msK: to model the effect of deeper, long-waveiengtn convection such as giant ceils that, we argue in 'jl.e, exchange 
very little energy with acoustic oscillations. 

Though, as Figure 1 -hows, it :s likely that the characteristic temporal and spatial scales of convection vary 
continuously across the convection zone, convective processes can be thought to be segregated into two concentric 
-neils fan outer shell and an inner sneil). with convection in each shell being dominated by distinct processes. In §I.e 
we attempt to quantify tne extent of each shell; here we discuss qualitatively the characteristics of the convection and 
:rs likely seismic effects in each -hell. 

The outer shell occupies rhe top few scale heights where the acoustic and convective physics are most complex. 
Convection is most vigorous in this shell, being highly turbulent and with relatively short characteristic convective 
lifetimes and length scales. In this shell, the convective velocity is an appreciable fraction of the local sound speed 
\M c: 0.3). the time scales of the turbulence and of the acoustic radiation are commensurable, and the amplitudes of 
the p-modes and the convective waves are largest. Goidreich and Keely (1977a, b) and Goldreich and Kumar (1933, 
1939) calculated the amplitudes and energies of the p-modes under the assumption they are excited by stochastic 
turbulent convection. Results of Goldreich and Kumar (1988) show that seismic wave emission and absorption in the 
Sun principally take place through interaction with turbulence in the top few scale heights of the convection zone. We 
ienne the radial extent of the outer shell as the region of significant interaction between the p-modes and convection. 
We argue in <jl.e that this region also marks t he extent of significant three-mode coupling and attempt to approximate 
:ts depth extent as well. 

In the outer shell, convective ceils evolve rapidly (Stein and Nordlund 1989; Title ti ai 1989). If, in addition, 
ceils are distributed isotropically in space, rhen they will produce little accumulated splitting effect on globally 
propagating waves. There vviii be local acoustic effects, but the isotropic assumption guarantees that the net global 
effect on frequency will be small. However, acoustic modal amplitudes, line-widths, and degenerate frequencies will be 
affected by outer shell processes [c. j. Brown 1984; Christensen-Dalsgaard and Fransden 1933: Christensen- Dalsgaard 
?t ai 1989) such as convectton-osciilation coupling, three-mode coupling (Kumar and Goldreich 1989), and radiative 
damping. 

The inner shell is much larger than the outer shell and lies directly beneath it. occupying, as we argue below, 
more than — 99.89c in radius of the convection zone. By definition, the emission and absorption of seismic waves 
by turbulence in this shell is negligible, and convection-oscillation coupling can be ignored accurately. Consequently, 
the anelastic condition can be applied. Furthermore, p-mode amplitudes are much smaller than in the outer shell 
and the solar gas in this shell is optically thick so radiative damping is negligible. Thus, the contribution to the 
interaction coefficient describing three-mode coupling in the inner shell is relatively small (Kumar and Goldreich 
1989). Therefore, we argue that splitting and the global distortion of seismic wavefronts dominantly result from 
convection that is relatively coherent temporally and spatially. If long-lived, long wavelength features of convection 
do exist, they would possess characteristic signatures in p-mode frequencies, line-widths, and displacement patterns 
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•hat could be computed from tlie theory presented herein. In principle, once we have identified these signatures, 
'heir observation wouid place constraints on the causative convective structures. For example, we show in a future 
cntnbution that because of the way heliosetsmic images are reduced and analyzed, the effect of aspherical structure 
s ro broaden line-widths. This effect --an G significant for modes with low intrinsic damping rates. The value oi 
:;eiioseismolo£ical constraints such as ihtse would be enhanced by the fact that large-scale convection has been linked 
with the dynamic structure of the differential rotation ( e.g r Gilman 1087) and with the solar dynamo (e.g. StLx 1981). 
!n audition, magnetic activity observed at the solar surface probably is controlled by us subsurface expression that is 
immateiy involved with Hows at depth. 

Next, we address two questions in uni I.e. (1) What is the evidence that large-scale convection exists in 

: h e inner shell? (2) What is the extent of t he outer shell where we do not accurately model the convective effect of 
convection? 


d) On the Existence of Giant- Cells 

A problem for the utility of the theory presented here is that giant cells have not been unambiguously observed 
u the soiar surface. If they do exist at the surface of the Sun. their amplitudes are less than 10 ms’ 1 (Howard and 
Labonte 1980: Labonte et ai 1981: Brown and Gilman 1984). Nevertheless, the evidence for their existence is strong, 
r housh circumstantial. (1) First, t he Sun displays a number of features that are suggestive of sustained large-scale 
motions [Gilman 1987). These include persistent large-scale patterns in the solar magnetic field, the coronal holes 
which survive several sohr rotation p^rious without being sheared apart by differential rotation, and the existence of 
active longitudes where new active regions preferentially arise. (2) Second, the observed distinct cellular convection 
:nav continue well below the surface. In the mixing length picture of convection (e g-. Fig. I). the scale of the 
'onvective eddies is set by tiie pressure scale-height so that one would expect layers of convection with monotonically 
increasing vertical scale with depth. In addition, both linear and nonlinear models (e.^.. Gough et ai 1976) have 
shown that even when the fluid is compressible, and the stratification includes several scale heights, convection 
spanning the enure unstable layer is favored. Thus, for the Sun, patterns of motion with horizontal dimensions up 
to the depth of the convection zone [i.e. A ~ 200.000 km or harmonic degrees of / ^ 20) wouid be expected. (3) 
Third, the space-lab experiment of thermal convection (Hart et ai 1986) and the numerical simulations of Glatzmaier 
1 1984) and Gilman and Miller (1986) have suggested that large and sustained patterns of motion may exist in the 
Sun with scales approaching the depth of the convection zone. (4) Fourth, Hill (1988) constructed three-dimensional 
spectra ( k z , k v , uj ) of helioseismic images of small rectangular regions near the solar equator and discovered relatively 
large-scale horizontal, poleward (lows of approximately 100 m/s that may be the surface expression of giant-cells. (5) 
Finally, a possible explanation of the small vertical velocities of the supergranules and the absence of a strong signature 
of giant-cells in the data of Howard and Labonte (1980) may be found in the work of Latour et ai (1981) and van 
Ballegooijen (1986). Latour e( ai ( 1981) found that buoyancy breaking in A-type stars may occur in upward-directed 
flows that have horizontal scales large compared to the pressure scale height of the region into which they penetrate. 
This leads to lateral deflection and strong horizontal shearing motions. If this result applies as well to G-type stars 
5 uch as the Sun. it may provide the explanation for the lack of surface observations of giant-cells. In addition, van 
Ballegooijen (1986) found that density stratification screens out periodic components of the near surtace flow pattern 
in his convection model so that periodic motions that exist at depth would not be observed at the surface. 

e) Justification of Linearization for Application to Giant Cell Convection 

We now attempt to quantify the extent of the outer shell, defined to be that region where convection-oscillation 
coupling is appreciable. The extent of energy exchange between oscillations and turbulent convection depends on their 
relative time and velocity scales. Perhaps the best available measure of the coupling between convection and acoustic 
oscillations is the flux of energy F r pumped into to the acoustic modes from the convective motions. Goldreich and 
Kumar (1989) derive an expression for F P given by 






( 12 ) 


where F c is the convective flux given by equation (11). F p depends primarily on the Reynolds stresses. 

An inspection of Figure 2. which plots the radial dependence of F pi reveals that convection-oscillation coupling is 
relatively insignificant below the top .15% of the convection zone. Thus, as a mechanism of oscillation-convection 
coupling, Reynolds stresses and entropy fluctuations act far more efficiently in the top few scale heights than in the 
deeper layers where the characteristic velocities are smaller and the length scales are larger. 
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Kumar and Goldreich (198'Jj also discuss rhe effect or' nonlinear interactions among solar acoustic modes. They 
-inue that these interactions are strongest in r he outermost layers of the Sun. Indeed, an inspection of their Figure 
:» indicates that the coupling coefficients are sensitive to three-mode interactions only in the outer — 2% by radius 
:{ *he convection zone. Consequently, we infer that rhree-mode interactions can be ignored in the determination of 
seismic erfecrs of convective Hows beiow this depth. 

In criticism ot the linearization in both convective velocity and seismic displacement, it misfit be suggested on 
intuitive grounds that a theory governing the eiFect of convection on seismic waves must be accurate along the entire 
path of the seismic wave, and since ail seismic waves propagate through the outer shell the theory must be general 
-nough to govern outer shell physics. This would certainly be true if we were interested in describing ail of the seismic 
erfects of convection. However, as discussed in turbulent convection and other nonlinear processes in the outer 
shell will dominantly affect modal amplitudes and degenerate frequencies. We are only interested here in determining 
split frequencies and perturbed eigenfunctions (or displacement patterns) of acoustic modes. Consequently, outer shell 
physics will be subsequently ignored. 

In conclusion, we define the outer shell to have a depth of ^ .2% of the convection zone and we argue that the 
seismic effect of convection can be modeled accurately with a linearized theory for flows within the the inner — 99.8% 
of the convection zone by radius. 


f) Overview 

In -jll we discuss reference frames and the separability of the equations of motion, and present a means of transfer- 
ring the theoretical results presented in this paper from a frame corotating with the Sun to an inertial frame that can 
be roughly identified as the observer s frame. In <>111 we present mathematical representations for convection and for 
the asphericities in the elastic-gravitational variables. The equations of motion governing the seismic oscillations in 
the presence of a steady-state global-scale velocity field and the associated static structural perturbations to density 
and bulk modulus are derived in *jIY. We derive in the quasi-degenerate perturbation theory needed to calculate 
the influence of a velocity held and structural perturbations on solar oscillations. In §VI, we derive the general ma- 
trix elements that determine the displacement field and split frequencies caused by an aneiastic model of convection 
represented with scalar and vector spherical harmonics by using the perturbatiou operator derived in §IV and the 
perturbation theory derived in $ V. The application of quasi-degenerate perturbation theory to the acoustic modes of 
the Sun and the derivation of the general matrix elements are presented in § VI. In §VH we discuss properties of the 
supermatrix. In § VIII. we consider differential rotation. In §IX we show how the the supermatrix may be used to 
generate theoretical seismograms. The principal conclusions of the paper are summarized in §X. 

The system of equations that governs the modal eigenfucntions and eigenfrequencies of the SNRNMAIS solar 
model is presented in Appendix A. The equation of motion of the non-SNRNMAIS solar model is derived in Appendix 
B. In Appendix C. we present a mathematical method adapted from Phinney and Burridge (1973) that considerably 
simplifies the application of differential operators to vector and tensor fields in a spherical geometry which are common 
in helioseismology. This technique is used to calculate the general matrices presented in §VI. Appendix D contains 
a discussion of the incorporation of the aneiastic condition into the general matrix and a derivation of the matrix 
elements for the Coriolis force, centripetal acceleration, and for general convection. Appendix E presents detailed 
expressions for the matrix elements for aspherical perturbations in the elastic-gravitational variables. 


[Sections II-IX omitted.] 

X. SUMMARY AND CONCLUSIONS 

The purpose of this paper has been to derive a theory that governs the effect of steady-state convection and 
associated asphericities in the elastic-gravitational variables (adiabatic bulk modulus k , density /?, and gravitational 
potential d) on seismic frequencies and displacement patterns and to present a formalism with which this theory can be 
applied computationally. The theory is not intended to predict nitMal amplitudes since these are governed, in part, by 
the exchange of energy between convection and seismic waves, which is excluded by our theory since it is linear and since 
the convective flow is defined to be aneiastic. To the best of our knowledge, every global-scale study of the heiioseismic 
effect of convection or structural asphericities, to date, has assumed an axisymmetric model. We have made no such 
assumption, and have represented convective flow (a vector field) and structural asphericities (a scalar field) with 
general global basis functions, vector and scalar spherical harmonics, respectively. We also represent the eigenfunctions 


ORIGINAL PAGE IS 
OF POOR QUALITY 




T rhe spherical reference model i Hie SNRXMAIS solar model) with vector spherical harmonics. These representations 
allow us to employ ouaM-dezeiieraift perturoanon theory in a straightforward manner to derive the general matrix 
•leinents that govern tin* modal coupling and splitting caused by convection and the structural asphericities. 

We present formulae for Hu? general matrix elements explicitly in terms of the scalar eigenfunctions of the SNRXMAIS 
-;iar model. Thus, the use of tins theory requires only the following quantities: { 1) a SNRXMAIS solar model U*(r) 
-iiu o\ r i > . >2) the -eismic scalar • ezenfunctions of the ^XRXMAIS solar model n 6’i(r), n \ H>i. n v'i(r), n doi{r) 

an a n^o.qrj). and (3) the spnericai harmonic representation of convection i if‘(r), u;(r). and U'j(r)) and/or asphericities 
:n the elastic-gravitational variables [dn : : [r) ami ep*(rj) at each radial knot of the SNRXMAIS model. The general 
matrix elements compose the llermitian supermatrix Z. whose eigenvalues are the eigenfrequencv perturbations of the 
general non-SNRNMAIS >olar model and whose eigenvector components are the expansion coefficients in the linear 
combination forming the perturbed eigenfunctions (or displacement patterns) in which SNRXMAIS eigenfunctions 
are basis functions. 

Optimally, the next stage of this research would be the application of this formalism to a realistic global model 
?f long- wavelength convection with associated asphericities in the elastic-gravitational variables. A major aspect of 
rais effort wouid be the determination of the accuracy of degenerate perturbation theory relative to quasi-degenerate 
perturbation theory. We have argued in 5 VII I that due to the spacing of modes that satisfy the selection rules, 
quasi-degenerate coupling r nn. for ail practical purposes, be neglected in modeling the effect of differential rotation 
m iielioseismic data. In -meet, modes that can couple through differential rotation are too far separated in frequency 
o coupie strongly. This i» not rhe case when we consider Honaxisymmetric flows and asphericities. In this case, near 
lezeneracies will regularly occur, and couplings can be relatively strong especially among SNRXMAIS modes within 
'lie same multiples However, since soiar convection is dominantly axisvmmecric, complete hybridization of modes 
would be rare, and the perturbed mooe would retain many of the characteristics ot a mode of the SNRXMAIS model. 
Most importantly, a perturbed eigenfunction, on average, wouid resemble a slightly perturbed SNRXMAIS eigenfunc- 
tion: t.e., a single spherical harmonic. However, these perturbations to the eigenfunctions and eigenfrequencies will 
be systematic and it should prove interesting to investigate the cumulative effect on splitting data. In particular, one 
could compute the perturbed eigenfrequencies for a given model of giant-cell convection (e.g., Gilman and Glatzmaier 
1981: Glatzmaier and Gilman 1931. 1932: Glatzmaier 1984) and then invert for the input differential rotation profile 
using currently standard methods that assume that nonaxisvmmetric components of flow are nonexistent. If the flow 
model were realistic, one would uncover any bias in the recently estimated differential rotation profiles. 

Another use of the theory would be to determine whether general asphericities in the elastic-gravitational variables 
could appreciably affect hclioseismic data. For example, Kuhn et al. (1988) observed a surface temperature variation 
of several degrees Celsius from the solar south pole to the solar equator and hypothesized that this or a similar 
structure may be responsible for the non-zero even-degree frequency splitting coefficients. Given an equation of state, 
these temperature variations could be expressed in terms of the perturbations Sp° 3 [r) and £«J(r). Although the depth 
extent of the observed temperature variation is unknown, different hypothesized depth structures could be constructed. 
Using the theory presented here, the general matrix elements and, hence, the splitting caused by each temperature 
model could be computed and used to test the hypothesis of Kuhn et al. (1988). 

The major constraint on the application of the theory presented here is that we have assumed that the convective 
flows and asphericities are stationary in time. Consequently, we view this paper as the first step toward a more 
general theory governing time-varying flows. Nevertheless, a number of the consequences of the theory will hold for 
time- varying flows as well. Most importantly, the selection rules listed in § VII will hold for non-stationary flows. 
For example, under self-coupling (or within degenerate perturbation theory), by Selection Rules 1 JC and 3 JC , only 
odd-degree s toroidal flows and even degree structural asphericities with s < 21 will affect the splitting and coupling 
of acoustic modes with harmonic degree /. 

In closing, since this paper is long, it is worthwhile to present a road map through the major results. Modal 
notation and terminology are discussed in §l.a and model notation and terminology are presented in §IV.b.i. The major 
assumptions of the theory are presented and discussed in §I.b, and are justified in and I.e. The mathematical 
representation of convection is in equations (23)-(25) and the representation of the elastic-gravitational variables is in 
equations (29)-(3I) and (35)-(37). The general equation of motiosi is equation (B . 1 1 ) and the equation of motion for 
the perturbed model with first-order perturbations including rotation, ellipticity in the structural variables, centripetal 
force, convective flow, and asphericities in the elastic-gravitational variables is equation (50). The general forms of 
the general matrix element and the supermatrix are shown in equations (67) and (68), respectively, and the general 
matrix element for the perturbations listed in the previous sentence is in equation (77). The explicit form of the 
general matrix elements suitable for computation, written in terms of the scalar eigenfunctions of the SNRXMAIS 


37 


ORIGINAL PAQ^ IS 

CF POOH QVALfTY 



^oiar model, can be found in equation (02) with notation and the integral kernels defined in equations \ U2). WV 

ronsider tills the mam result ot r h is paper. Three selection ruies govern in 2 ; coupling are listed in equations ( 120). ( 122) 
. and 1 122). with the seif-coupling term 01 the selection rules inequations 1 1 J L ? . (122). and f 124). The Diagonal Sum 
Rule and the huperaiagonal Sum Rule are star.-i and proved in 'jYII.d. The general matrix element and selection rules 
:or differential rotation are in equations 1 127} unci ■ 141), respectively. All results of the paper are presented in a frame 
corotating with the Sun. Equation i Is) can be used to construct the perturoed eigenfrequencies and eigenfunctions 
ci a non-SNRNMAIS solar model in an inertial frame. The Generalized Spherical Harmonic formalism, which was 
used to derive the explicit form ot the general matrix elements, is discussed in Appendix C. The incorporation of 
inelasticity constraint into the general matrix element for convection is the subject of Appendix D. 
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Figure. 1. - Characteristic length, velocity, and time scales of convective eddies a 5 predicted by mixing length 
theory (eqs. [7], [9] and [10]) plotted as a function of depth. The solar model of Podsiadlowski (1939) was used to 
calculate these quantities. 
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Figure. 2, - A plot of F p M j F e 'ahich is the flux of energy pumped into the acoustic modes from the convective 
motions (see eq. [12]). We have normalized F p by its peak value. We use the radial dependence of F P to argue that 
coupling between convection and p modes is significant only in the top ~ 0.2% of the convection zone. 
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ABSTRACT 

We consider wave generation by turbulent convection in a plane parallel, stratified atmosphere that sits in a 
gravitational field, g. The atmosphere consists of two semi-infinite layers, the lower adiabatic and polytropic 
and the upper isothermal. The adiabatic layer supports a convective energy flux given by mixing length 
theory; F c pv \ , where p is mass density and v H is the velocity of the energy bearing turbulent eddies. 

Acoustic waves with w > a ) ac and gravity waves with co < 2 k h Hi(o b propagate in the isothermal layer whose 
acoustic cutoff frequency, ov, and Brunt- Vaisala frequency, c o ft , satisfy co ^ = yg/4Hi and col = (y — \)g/yH h 
where y and H i denote the adiabatic index and scale height. The atmosphere traps acoustic waves in upper 
part of the adiabatic layer (p-modes) and gravity waves on the interface between the adiabatic and isothermal 
layers (/-modes). These modes obey the dispersion relation 



for co < a)". Here, m is the poly tropic index, k h is the magnitude of the horizontal wave vector, and n is the 
number of nodes in the radial displacement eigenfunction ; n — 0 for /-modes. 

Wave generation is concentrated at the top of the convection zone since the turbulent Mach number, M = 
v H /c, peaks there; we assume M t < 1. The dimensionless efficiency, rj, for the conversion of the energy carried 
by convection into wave energy is calculated to be rj ~ Af/ 5/2 for p-modes, /-modes, and propagating acoustic 
waves, and r\ — M t for propagating gravity waves. Most of the energy going into p-modes, /-modes, and prop- 
agating acoustic waves is emitted by inertial range eddies of size h ~~ A f r 3/2 H ( at co ~ and k h ** 1 /H r The 
energy emission into propagating gravity waves is dominated by energy bearing eddies of size and is 
concentrated at co ~~ vJH, ~~ M t co,*. and k h ~ l/ft v 

We find the power input to individual p-modes, to vary as tt / 2 » J + 7 "- 3 >/<" +3 > at frequencies co < vJH v 
Libbrecht has shown that the amplitudes and linewidths of the solar p-modes imply £ p cc co 8 for 
co 2 x 10” 2 s“ l . The theoretical exponent matches the observational one for m % 4, a value obtained from 
the density profile in the upper part of the solar convection zone. This agreement supports the hypothesis that 
the solar p-modes are stochastically excited by turbulent convection. 

Subject headings: convection — Sun: atmosphere — Sun: oscillations — turbulence — wave motions 

I. INTRODUCTION 

Lighthill (1952) wrote the seminal paper on the generation of acoustic waves by turbulence in homogeneous fluids. Stein (1967) 
extended LighthilFs techniques to stratified fluids and also treated the emission of gravity waves. We reconsider Stein’s problem for 
a more realistic model atmosphere and relate the turbulent spectrum to the convective energy flux via the Kolmogorov scaling and 
the mixing length hypothesis. Our goal is to estimate efficiencies for the conversion of the convective energy flux into both trapped 
and propagating waves. Wc treat mode excitation but not mode damping. Thus, we cannot estimate the energies of trapped modes 
which depend upon the balance between these two effects. 

The plan of our paper is as follows. In § II we describe the model atmosphere and its eigenmodes. Next, in § III, we derive 
expressions for the rates at which individual modes gain energy from turbulent convection. In § IV, we estimate the total emissivities 
for the different wave types, p-modes, /-modes, propagating acoustic waves, and propagating gravity waves. A comparison of our 
results with those obtained in earlier studies, and a discussion of their implications, is given in § V. 

n. ATMOSPHERE AND EIGENMODES 
a) Static Atmosphere 

Our model atmosphere is plane parallel, sits in a constant gravitational field, g, and consists of two semi-infinite layers, the lower 
adiabatic and polytropic and the upper isothermal. The pressure, p, density, p, and temperature, T, are continuous across the 
interface between the two layers. In the lower layer the adiabatic and polytropic indices are related by f = 1 + 1/m. The adiabatic 
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index in the upper layer, y, may differ from f\ The z coordinate measures depth below the level at which the adiabatic layer would 
terminate in the absence of the isothermal layer. We denote quantities evaluated at the top of the adiabatic layer by a subscript t. 
Parameters in the isothermal layer are distinguished by a subscript i. Note, the ratio of the sound speeds cjc t = (y/F) 1/2 . 

In the adiabatic layer the thermodynamic variables exhibit a power-law behavior with depth: 



The sound speed, c, and the pressure scale height, //, satisfy c 2 — gz/m and H = z/(m 4- !). 

The isothermal atmosphere is still simpler: T = T t , c = c i9 and H = H ( are all constant, whereas p and p are proportional to 
exp {z/ Hi). 

b) Normal Modes 

We choose the Eulerian enthalpy perturbation, Q = p x /p y as the dependent variable in the linear wave equations. These read 


in the adiabatic layer, and 



( 2 ) 

( 3 ) 


in the isothermal layer (Kumar and Goldreich 1989). Here, co is the wave frequency and k h is the horizontal wavevector (/ = k h R Q ). 
The displacement vector, <*, is related to Q by 


in the adiabatic layer, and by 





Si 


J_£Q 

co 1 cz 



1 p Q } (y- 1) 

(a ) 2 -<ot)_dz yHi 



( 4 ) 


( 5 ) 


in the isothermal layer. 

The normal modes are obtained by solving equations (2) and (3) subject to Q -► 0 as z -► oo, Q and continuous across the 
interface at z f , and the appropriate boundary conditions as z -► — oo. The continuity of £ k follows from that of Q. 

The modes are classified as trapped or propagating, and as composed of acoustic or gravity waves. The adiabatic layer supports 
acoustic waves, but not gravity waves. Moreover, it refracts acoustic waves upward. Thus, propagating modes must be traveling 
waves in the isothermal atmosphere. 

Solutions of the wave equation in the isothermal atmosphere are proportional to exp ( — k ± z), with 





where co M and co b are the acoustic cutoff and Brunt-Vaisaia frequencies: 


( 6 ) 


(U 


2 

ac 



( 7 ) 


and 


2 (y — l )g 

cot =- 




( 8 ) 


Thus ml = My - i)o£/y 2 . There are two branches to the dispersion curve for traveling waves. For 2k h H { 4 1, these are a high 
frequency, acoustic wave, branch with co > and a low-frequency, gravity-wave, branch with co < 2k h H { c o b . 

Wave excitation by turbulent convection is concentrated in the upper adiabatic layer where the convective velocity peaks. We 
seek analytic expressions for the normalized eigenfunctions in this region. Since the dominant interactions are proportional to 
d 2 Q/dz 2 (see § Illh), we explicitly evaluate this quantity for each mode. In doing so, we drop factors of order unity including, in 
places, y, T, and m . 


i) Trapped Modes 


Trapped modes correspond to evanescent solutions in the isothermal layer and are restricted to a discrete set of eigenfrequencies 
for fixed k h . In the limit that the adiabatic layer extends to vanishing surface pressure, the eigenfunctions may be expressed in terms 
of associated Laguerre polynomials and the dispersion relation reads 



( 9 ) 


3 . 
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Christensen-Dalsgaard and Gough 1980)^rapp^ eigenfunction (Christensen-Dalsgaard 1980; 

gravity waves, or /-modes. Trapped o-modes with n * () Hn nnf l-, ?- u cor ’’? s P ond to n 0. Modes with n = 0 are surface 
Brunt-Vaisala frequency vanishes. Equation (9) remains^ eood annroTtW 06 ad,aba,lc ]a y er is neutrally stratified, that is, its 
Only the physical solution, the one that u < “V eVen With fin,te surface P ressure - 

condition reads g ° WS l£SS rap,dly Wlth hei « ht In the >sothermai layer, is normalizable. The normalization 


1 s I • «• = 6 m . , 

J - oo 


. - ( 10 ) 


* dO 

- dz 4 
Jt, c 


Q»QZ = K* 


(ii) 


. I. P- Modes 

It is easily shown that 


and 


2kV 

co J 




(12) 


(13) 


equation (2) dominates and Q oc exp £Vanescent ' Both Q and $ increase slowly with height above z, . Below z 2 the k h term in 

We study the p-mode eigenfunctions in the dual limit o « m and Jh « ^ i i 
pressure the eigenfunctions are solutions of equation (2) that are"analvtic a/t' -nl P ° 1 1 ytrop,c layer with vanishing surface 
assoaated Laguerre polynomials. When the polytropic layer is ov«W kJ 7 ° : S ° 1 u , tIons may 1x5 expressed in terms of 

contribution l from the solution that is singular at z = 0. However the Sll , ,s .° hcrmal la y er - the eigenfunctions include a 
ensure that the contribution from the singular solution is small for o> « m ^ on 1 ,on s at the interface between the two layers 
We can approximate the eigenfunction in the region of propagatfon^z^^ z < z 2 , by the WKB solution 


Q~ P 


("i- l)/2 


B p sm 


iin [Ht)"’ + *>] • 


(14) 


zone above .he atmosphere 

» ^oTors,K — - — *«. - 

i tsszr. s s ~ 

below Z . z,. However, this is of httle consequence for the eigenfunctions that become evaneSent well 

Given the properties of the eigenfunction desert hed above, it follows from the normaliaation equation (1 1) that 


b;~; 


w2(m - l 




Evaluating d 2 Q/dz 2 we obtain 


forz r < z z x . 


g {m ~ 2) p t 




(15) 


(16) 


, 2. f- Modes 

equations in both the adiabaficand^ ? P fonmld f W '‘ h ^ = ^ <S a " exact solution of ,he w *ve 

Thu family of normal modes consists of gravity waves CMfit^^ear the' smthtto^the conwwrton zone? theV are known^s/modes! 
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The /-modes are incompressible, V * f = 0, which accounts for their simple dispersion relation. The amplitude, B f , is determined 
from the normalization equation (1 1) to be 


For ail r. 


B/ 


:?a> 2lm - l) k„ 


(17) 


c 2 Q 


= B f exp { — k h z) . 


( 18 ) 


ii) Propagating Waves 

Modes that propagate in the isothermal layer have continuous spectra. They are chosen to have no net flux in the isothermal 
layer; that is, they are composed of pairs of inward- and outward-propagating waves of equal amplitude. This choice ensures that 
propagating modes have real frequencies and are orthogonal to trapped modes. These modes are normalized such that 


dzpi„ * = S(to- to ') , 


( 19 ) 


at fixed k h . The upper limit on the integral in equation (19) may be taken to be r |t since the contribution from the adiabatic layer is 
finite, and therefore negligible. 

!. Acoustic Waves 

These modes have to > and propagate in the isothermal atmosphere and in the upper part of the adiabatic layer. They are 
evanescent below the lower turning point at z 2 ~~ to 2 /gk b . We deduce the properties of the eigenfunctions in the joint limit to P to x 
and k h <| a>/c r *. 

In the isothermal layer 


Q = C a sin [KAz, - z) + £ J exp > 

where K z ^ cu/c r . Application of the normalization condition given by equation (19) to equation (20) yields 


( 20 ) 


C 2 

^ a 




P I 


( 21 ) 


We approximate the eigenfunctions in the adiabatic layer by the WKB solutions 


Q 




( 22 ) 


for < z « z 2 . The continuity of Q and £. across z, is used to relate B„ and 4>. to C e and £„. The phase, 4>„ is determined by the 
condition that Q oc exp ( — k h z) for z -* oo. For a> just above B a (u), k h ) displays sharp ridges along extensions of the p-mode 
dispersion curves. These correspond to resonances for the scattering of incoming waves by the atmosphere. These ridges flatten for 

uj 5> (b k and 




TC] 

[1 + (r - y) cos 2 4>J 


For later use we record 


(23) 


c 2 Q 

dz 2 


mco 2 „ 

— C, 


(24) 


for z, < z <$ z 2 . 

2. Gravity Waves 

Gravity modes with to < 2k fl H i (o b propagate in the isothermal atmosphere but are evanescent in the neutrally stable adiabatic 
layer. We detail their properties in the double limit to <$ 2k b //, to b and 2k h Hi < 1. 

In the isothermal atmosphere 


Q = 


C g sin [K r (z, - z) + CJ exp 



where K z ^ i(o h /to)k h . 

The amplitude C 9 is determined by the normalization equation (19) to be 


C 


2 

9 



(25) 


(26) 
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miv r^'„ f ^H 2 ; f ‘ * k V thel u term 10 the wave equatlon (2) is much smaller than the first and second tenns and 

may oe ignored, i he reduced wave equation yields 


7 


■r \(W- n 


+ D a 


(27) 


The ratio D g ,B g is determined by fitting Q in the upper part of the adiabatic layer to Q x exp ( — k h z) at r k l . For 2Jt // <^ 1 

g/ &g ^ 1 • H t 1 

The continuity of Q and c z across z, is used to relate B g and D g to C and ( . We find 


9^9 g' 

„ o)k h zf' 2 
tan t. 177^ • 


(28) 


The small value of tan ( is due to the change in orientation of the velocity field from almost vertical in the top of the adiabatic laver 
to almost horizontal in the isothermal layer. From equation (28) it follows that y 


s, 2 


z^co% 

9 1!2 Pi 


For later use we note that 


holds for z <1 k h l . 


c 2 Q njm - \)Q 


(29) 


(30) 


c ) t urbulent Convection 

In the absence of a reliable theory for turbulent convection, we are guided by the mixing length hypothesis. According to this 
Sht ,S f carned fl by turbulent wb °se dimensions are of order the local pressure scale 

related t^ the"mean e^tropy^adicnt! ds/dz! |^y tr °^ y flUCtUat, ° nS aSSOC,3,ed ^ th ~ b “ ri “« edd “* «**> a " d ^ a " 


gjp ds 
c„ dz 


where c p is the specific heat at constant pressure per unit mass, and 


These relations lead to 


"s- 


(31) 


(32) 


Since F c is independent of z. 


F C - pTv„s„ ~ pv 3 „ . 


(33) 



(34) 


where v, = r„(z,). 

In treating the convection zone as adiabatic we have been neglecting the superadiabat.city of the 
L " Tds/dz, with respect to the adiabatic temperature gradient, g/c„. From equation (32) it follows that the 

m \j hp PvnraccAH up r ’ 


temperature gradient, 
ratio of these gradients 


T ds 
g dz 


M 2 , 


where the Mach number of the turbulence, M = v H /c. Appeal to equation (33) establishes that 


(35) 


M 



(36) 


^ aS * U ™ c ‘ hat ,hc turbule " t velocities are substantially subsonic even near the top of the convection zone, that is, M, « 1 . Under 
modes dltl ° nS WC 3re JUSUfied in a PP roximat,n g the convection zone as adiabatic when calculating eigenfunctions for the normal 

The characteristic time scale of the energy bearing eddies is 




V H 


(37) 
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It is smallest at the top of the convection zone where 


1 

M t “>ac 


(38) 


The velocities of smaller, h < H< inertial range eddies are related to those of the energy bearing eddies by the Kolmogorov scaling 
(Tennekesand Lumley 1972), 



(39) 


at fixed c. The Kolmogorov spectrum applies to turbulent convection because, below the scale of the energy bearing eddies, the 
Reynolds stress provides greater accelerations than the buoyancy forces (Goldreich and Keeley 1977a). This implies that entropy 
mixes like a passive scalar contaminant in the inertial range. Thus, 


S H 



(40) 


The depth dependence of the properties of eddies of fixed size h follows from equations (32), (34), (37), and (40). We find 

r hz? i *' 3 

~ y rl ^Tm)J 

r z {m+i) h 2 i i 
x ^) - T { ] 

S*(z) 


11/3 


(41) 


hz\ 2m + 3> 1/3 

S f ,2<m + 2) 


III. MODE EXCITATION 
a) Source Terms 

We begin this section by adding source terms due to turbulent convection to the linear wave equation (2) for the adiabatic layer. 
Next, we classify the individual terms as sources of monopole, dipole, and quadrupole radiation. Then we evaluate the excitation of 
wave modes by these sources. 

We distinguish three principal sources of wave excitation by turbulent convection. They are, the expansion and contraction of 
fluid due to the gain and loss of specific entropy, buoyancy force variations associated with these entropy changes, and momentum 
transport by the fluctuating Reynold’s stress. 

We derive the inhomogeneous wave equation from the linearized versions of the equations for mass and momentum conservation 
supplemented by the equation of state for a perfect adiabatic gas. We augment the momentum equation by the divergence of the 
turbulent Reynolds stress, and the adiabatic equation of state by the entropy fluctuations associated with turbulent convection. 
These equations now read : 

+ V •(/>») = 0 , (42) 

ct 

— + Vp, - p x g = - V • (pvv) = F , (43) 

Ct 


Pi __ £pi = ± 

P P ?v' 


(44) 


where /?,, r, and s are the Eulerian density, pressure, velocity, and entropy perturbations associated with the turbulent 
convection and the waves it generates. The subscript 1 attached to the density and pressure perturbations denotes that only the 
lowest order variations of these quantities need be retained. Equation (44), the Eulerian form of the perturbed equation of state, 
holds because the background state is isentropic. 

Eliminating p l and v from the left-hand sides of equations (42H44), we obtain the inhomogeneous wave equation 


££_!<.!£_ s <l> + s<2) 

V c 1 dz c 1 dt 2 p 


(45) ' 


where 



(46) 


6 


S' 2) = V ■ F . 



700 


GOLDREICH AND KUMAR 


Vol. 363 


The interpretation of equation (45) is somewhat subtle. Provided we drop the final c~ 1 d 1 Qlct 1 term on the left-hand side as a first 
thC ° fsubso “f tur . b , u if"“’ '* determines the near field turbulent pressure perturbations from the turbulent 
i y d cntropy perturbations. The c 5 Q/dt term connects the near field perturbations to the wave field perturbations The 
latter may be expanded in terms of the normal modes. 

The identification of sources by multipole order is a useful device in estimating wave emission by turbulent convection It helps to 
separate the sources that must be retained from those that may be safely discarded. For homogeneous and isotropic turbulence the 
multipole expansion may be earned out in several equivalent ways. In our application the turbulence is r-dependent, and therefore 
inhomogeneous, and the atmosphere is stratified, and therefore anisotropic. Under these circumstances the method of choice is to 
identify sources according to whether they involve a change in fluid volume (monopole terms), a source of external momentum 
(dipole terms) or merely internal stresses (quadrupole terms). 2 Classification based on the angular dependence of the wave 
amphtude m the radiation zone is not useful, because the angular dependence results, in part, from the anisotropy of the medium. 3 4 

Identification of sources by the number of their spatial derivatives also leads to ambiguity, since it differs according to the choice of 
dependent variable. 

ThI h ^ St ^ nn m ^ 1 )™? dire ? ly u from the v , oiume change due to an entropy change at fixed pressure. It is a monopole source. 
The second termini'' reflects the buoyancy force variation associated with this volume change. It involves a variation of the 
density of momentum supplied by the external gravitational force and is a dipole source. The double divergence of the Reynolds 
stress in S' reflects the redistribution of momentum by internal stresses. It is a quadrupole source. 

One might suspect that the monopole and dipole terms in S' 11 produce more acoustic radiation than the quadrupole term in S' 2 ' 
Treating these three terms independently appears to confirm this suspicion; the monopole and dipole terms are found to excite 
comparably greater amounts of acoustic radiation than the quadrupole term. However, the correct solution is more subtle As we 
qSupd^emSon dCStrUCtlVC ,nterfercnce causes thc total monopole plus dipole acoustic emission to be of the same order as the 


b) Amplitude Equation 

Hie total enthalpy perturbation, Q(x , r), is expanded in terms of the normal modes, Q a (z ), as 


Q = ~ ?= X lA.Q* «P ( ~io)t + ik h -x) + A • Q; exp (iajt - ik„ ■ x)] 


(47) 


7, h rfJJ S *, h ! CTOSS k seCtion °f thc atmosphere." The mode amplitudes, AM are slowly varying functions of time, 

\dAJdt\ « a\A,\. Substituting this expansion into equation (45), multiplying both sides by Q* exp (ion - ik h • x), and integrating 
over space and time, we obtain k 8 8 


AJi t) = 


i f r 

2icuj* rl;2 _ J 


d i xQ*(S tu + S' 21 ) exp (icor — ik h • x) . 


(48) 


Taking oo for the lower limit on the integral over t involves the implicit assumption that damping erases the memory of 
excitations from the distant past. 3 

Next, we integrate by parts to transfer all time and space derivatives to the eigenfunctions. The contributions due to the individual 
source terms are discussed separately below. 

The monopole plus dipole terms contribute 


^'(t) = 


1 


dt 


d i y 


ps 


K-f) 


2ia)s/ lli 

With the aid of the homogeneous wave equation (2), we transform equation (49) to 


exp (ion - ik h • x) 


^ * " 2iW 12 J..* j di X 17 Q*) exp {icot ~ ik " ‘ • 


The contribution due to the quadrupole term is 

A[ 2 \t) % 


2 iajjtf Lil 


dt 


d i xpvv:V\Q* exp (ion - ik » • x) . 


(49) 


(50) 


(51) 


The normal mode eigenfunctions share the property that kj Q. \ < \cQJdz\ near the top of the adiabatic layer. More precisely, 
other than the /-modes for which dQJdz = - k h Q„ the approximate mode eigenfunctions calculated in § life satisfy the strict 
inequality. This implies that 


A"\t) 


I 


liajs/ 


mf *1 

J - n m 


pc 2 s d 2 Q * 

d x T-y exp {icot - ik h • x) , 

c p cz 


( 52 ) 


2 This method preserves the ordering of source terms by the efficiency with which they generate radiation. 

3 For example, a spherically symmetric point source radiates amsotropically in a stratified atmosphere. 

4 For the moment we are treating the atmosphere as being of finite horizontal extent. 
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and 




dt 


d 3 xpv: 


2 c 2 q: 


cz 


exp (icot - ik H ■ x ) , 


(53) 


provide order of magnitude estimates for ,4', n (r) and A {2) (t). However, = 0 for /-modes as a consequence of their incompress- 
ibility. ...... 

Now we compare the relative sizes of A[ l) (t) and A i2 \t), We start with the contributions made by energy bearing eddies and go on 

to investigate those due to smaller, inertial range eddies. 

According to equations (31H32), c 2 s„/c p - v 2 „. Thus, except for the /-modes, the entropy and the Reynolds stress sources 
associated with energy bearing eddies make comparable contributions to A 2 (t ). This illustrates the destructive interference between 
the monopole and dipole amplitudes to which we referred earlier; for energy bearing eddies and acoustic modes with <o - v a /H, the 
monopole and dipole terms in equation (49) are each larger by a factor ^(c/u w ) 2 than the combined term in equation (50). The 
destructive interference between monopole and dipole amplitudes is a consequence of the anisotropy of the adiabatic layer. This is 
expressed by the anisotropic form of equation (2) which transforms equation (49) into equation (50). 

For inertial range eddies, c 2 sjc p v^(H jh) 1 ^ . This suggests that, unlike energy bearing eddies, inertial range eddies might excite 

waves more by their entropy sources than by their Reynolds stress sources. In fact, this is not the case. F rom equation (50) we see 
that wave excitation by the entropy source depends upon the time variability of the Eulerian entropy field. Inertial range eddies 
contribute to this time variation in different ways. The kinetic energy in an eddy of size h < H may dissipate raising the local value of 
s h . Neighboring eddies of similar size having opposite signs of s h may collide and mix their fluid thereby smoothing the spatial 
variation of the entropy field on scale h . An eddy of size h carrying an entropy fluctuation s* may be advected at speeds up to v B . Of 
these possibilities, the dissipation of kinetic energy into heat produces the largest entropy source. However, this source is just equal 
to that provided by the Reynolds stress. Thus, from here on we use equation (53) to estimate the total excitation rate of normal 

modes. .... . 

Destructive interference between monopole and dipole radiation fields holds the acoustic emissivity of turbulent convection at the 
level characteristic of free turbulence 5 for which the emissivity is dominated by acoustic quadruples. We did not appreciate this 
point in our earlier treatment of acoustic emission by turbulent fluids (Goldreich and Kumar 1988). There we discussed the 
emissivity of turbulent pseudo-convection, a surrogate for turbulent convection. Since this model has acoustic dipoles but not 
acoustic monopoles, its emissivity is greater than that of free turbulence. 


c) Excitation Rate 


Turbulent convection consists of a hierarchy of critically damped eddies. Different eddies of similar size are assumed to be 
uncorrelated. This assumption enables us to divide into several steps the calculation of the rate at which turbulent convection 
pumps energy into a wave mode. 

To begin, we estimate the magnitude of the incremental amplitude, A A\ y produced by a single eddy of size h located at depth z 
over its lifetime h/v h . 


A A\ - 


pv H h* c 2 Q* 
lioijrf 112 cz 2 


CO < T* 1 


(54) 


In arriving at the above equation we have assumed that the eigenfunction does not vary dramatically over Az = h ^ H. This is a 
good approximation for all the modes we are concerned with. At frequencies much greater than t* \ A A* declines exponentially 
with increasing cu. 

Next, we note that 




1^:1 


(55) 


is the mean rate at which one eddy supplies energy to mode x 
Then, summing over eddies of all sizes and depths, we obtain 


a> 


dzp 4 


r?**- 


(56) 


where 


(-*) 


Hiz) 

1 + [WT„(2)] 3/2 ‘ 


(57) 


In deriving equation (56) from (55), we include a factor sJdz/h\ the number of eddies in the horizontal slice of cross-sectional area sf 
between vertical depths : and : + dz. The appearance of dh/h in equation (56) denotes that each inertial range eddy accounts for a 
finite range of scale size dh/h ~ 1 . Carrying out the integration over h yields 






(58) 


? Free turbulence is turbulence that is not subject to external forces. 
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where the weight factor, W , is given by 


The weight factor is sharply peaked about 


where 


mu) = 


[1 + {toT r ) 3 ' 2 u ( 


15 * 


U* ^ I + 


(0>T f ) 


,3/<« 


(59) 


(60) 


* (WT,) 3( "’ +4) ' , ”'" 3 »[l + (Ci)T r ) 3 * 3m + 3 * m + 3 *J 1 (61) 

it decays as u lm + A) for Mw, and as u“ (3m ^ 7)/2 for u > u,. 

The peak in W is so sharp that £, is dominated by contributions from z ~~ for all wave modes. Physically, this means that the 
excitation is concentrated in the layer where the turnover time of the energy bearing eddies is most nearly equal to the mode period. 
This enables us to further simplify the expression for £, to 

£ ^P±Hl 1 

1 X, (CJT,)* 3l " * 3 3 *[1 4 " (CUT ( ) 3 * 3 — + 3 ^ 3 *— + 3 *] 




dz 2 


(62) 


IV. FLUXES OF ENERGY 

To evaluate the total excitation rate for each type of mode, we substitute the relevant expression for c 2 Q a (z*)/dz 2 given in § 116) 
into equation (62). Following that, we integrate £, over all modes of the family to determine the fraction of the convective energy 
flux that family receives. 

The frequencies of trapped modes satisfy equation (9). The flux of energy going into modes of a given family is 

(63) 

where the sum over a includes all modes in the family, the sum over n includes all dispersion ridges in the family, and j dk H is over all 
modes along a_ridge. The last equality follows because the spacing between adjacent k h modes in a box of horizontal area, s4 y is 
equal to Injjrf. Therefore, the number of modes in d 2 k h is stfd 2 kj(2n) 2 = (stf/2n)dk h k H . 

For propagating modes, to and k h are independently specified. The flux of energy into a family of modes is computed from 

F t -±Z£ = ±jdo>jdk h k h £', (64) 

where the double integral is over all modes in the family. 


a) P- Modes 


From equations (16) and (62), we obtain 


£ p ~ P, H? I’i 3 7-7-7 


+ *7 m - 3>/<« + 3) 


(65) 


1 + (a>T,) 3 ‘ 3 "' f5>/J, "- J ' ' 

At fixed k h , E p varies as + for cox, < 1 and as to 14 " - "’’' 2 for a>x, > 1. Toobtain the energy input rate per mode along 

the n'th p-mode ridge, we eliminate k h from equation (65) by using equation (9). This procedure yields 

* p, Hf f 3 Af? ,m * 31 (ore,)* 3 ’" 2 + 3>/l " + _ 3> 

p ~ (n + m/2) 1 + 5 » /3 «" +3 ' ‘ 

The total flux of energy going into p-modes follows from substituting equation (66) into equation (63): 

F, - 

From equation (66), we note that for tor, £> T the energy input rafeTs proportional to(tor f )* 4 " ,_3>/2 , which increases with increasing to 
for m > Since the maximum frequency for trapped p-modes is to^ , most of the energy flux goes into modes whose frequencies lie 
just below the acoustic cutoflf, to < w ac , and is emitted by inertial range eddies with h ~~ Mf n H t located in the top scale height of the 
convection zone. 

_ _ b) F- Modes 

The calculations for the /-modes are similar to those for the p-modes. We substitute equation (18) into equation (62) and find 

+ 7m- 3>/<m+3) 


( 66 ) 


(67) 


(C0T f f 


r ^ n H 3 r 3 Xf 2{m + 2) k 

iVi t ** I + (0)t f ) 3|3 '" + 5>/2| '" f3) 


( 68 ) 
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where we have set exp (-k h :J ~ 1 since (c,;:,,, ~ iVf, 2 (coT I j l2 " +3l/< "”' 3 '[l + (a>T I ) 3, " , + 3) ] < 1 forco < u > The rate of energy input per 
mode along the /-mode ndge reads 




2(wi + 3) 


(COT ) (2,w2 + 9m + 3V(in*t* 3) 

1 +(wt i ) 3<3 " , " 5, ' 2, ’ ,+3 ’ ' 


(69) 


The total flux of energy going into /-modes is 


Ft - P,v?M } s ' 2 = M ( 15 , 2 F C 


(70) 


From equations (24) and (62), we obtain 


c) Acoustic Waves 


£ a ~p,n?vf 


M 3 


(WT,)‘ 


(71) 


after averaging over the phase <t>„. Substituting equation (71) into equation (64), we derive the total flux of energy carried by the 
acoustic waves: 


F m ~ p,v}M } 512 = M} il 2 F c . 


(72) 


Most of this energy is emitted by inertial eddies of size h < Xf? l2 H ; located in the top scale height of the convection zone. It is carried 
by waves with to > to^ and k H < 1 jH r 


Equations (30) and (62) yield 


d) Gravity Waves 


£,~p,HtvfM t k h 




~ 3 ){(m + 3) 


[1 + (COT,) 


3(7« + 9)/2(m + 3 


r 


(73) 


so the power input into gravity waves peaks for cui f ^ 1. Equation (73) holds for k h in the range co/co b < 2 k b H t < (coT,) 3/(m+3) / 
[1 4- (cuT f ) 3y<mi ’ 3 *j. Substituting equation (73) into equation (64), we find the total flux of energy carried by the gravity waves: 

F 9 ~~p t v?M t = M t F c . (74) 

Most of this energy is emitted by energy bearing eddies located in the top scale height of the convection zone. It is carried by waves 
with cut, ~ l and k h < \/H r The vertical wave vector of these waves in the isothermal layer is k z \/(M t H ( ). 


v. DISCUSSION 
a) Previous Results 

Our principal results are dimensional efficiencies, r\, for the conversion of the convective energy flux into the energy flux in 
different types of wave modes; r\ f ~ r\ a ~~ M f l5/2 , and r\ q - M r It is illuminating to compare these efficiencies to those obtained 

in previous investigations. 

The classic result for the efficiency of emission of acoustic waves by homogeneous, isotropic turbulence is that of Lighthill (1952). 
Translated into our notation it is t] a - Af r 5 . Here we are thinking of the acoustic emission from a layer of turbulent fluid of thickness, 
H n embedded in an otherwise uniform atmosphere. The energy bearing eddies are characterized by size, H n and velocity, v r In this 
system, the acoustic emission is dominated by the energy bearing eddies, and is concentrated at a > ^ vjH t% k ~ M t H r We find 
*1 P ~~ n a ^ M, I5/2 , with the emission dominated by inertial range eddies of size h Mf l2 H t and concentrated atcu ~ cjH t , k h ~~ 1 /H v 
There are two relevant comparisons between our results and those of Lighthill. 

First, we can redo the estimate for from Lighthill's treatment restricting attention to emission from inertial range eddies having 
h < M? n H r These eddies, whose lifetimes x h < to’ 1 , dominate the emission of energy into p-modes and acoustic waves in the 
stratified atmosphere. A simple calculation yields rj a - M} 5i2 . This result agrees with ours showing that the acoustic emission from 
eddies with h < M? /2 H t is not affected by stratification. 

Second, we can modify our calculation of r\ p so that only the emission by energy bearing eddies is included. This is accomplished 
by repeating the procedure described in § IVu) but now limiting the integration over frequency along the p-mode ridges to 
a) < vJH r This exercise yields r\ p - M}°. The factor Mf by which this result differs from Lighthill’s may be accounted for as follows. 
Both in a homogeneous atmosphere and in our stratified atmosphere, the acoustic emissivity is proportional to | WQ | 2 . However, 
for to ~~ v t /H n | VVQ | 2 — (MJH t ) A \ Q | 2 in the homogeneous atmosphere, whereas | W Q | 2 (M?/H t ) 4 1 Q \ 2 in the stratified atmo- 
sphere. This difference, which accounts for four factors of M f , arises because p-modes with co ~ vJH t ~ are evanescent near 

the top of the convection zone in the stratified atmosphere, 6 The fifth factor of M t arises from differences in phase space mode 
densities. In a uniform atmosphere, the number density of modes having co v t /H ( is approximately (jVf f /H r ) 3 . This becomes Mf/Hf 
per unit area for a layer H t thick. The corresponding area density of p-modes in the stratified atmosphere is Mf/H 2 , just one power 
of Af f smaller. 


For acoustic waves with oj Z (o* . 1 VVQ | 2 is of the same order in the stratified atmosphere as in a homogeneous atmosphere. 
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